# Heavy vs. light flowers
# Callin Switzer
# 12 October 2016
# update 8 Dec 2017
# 
# Compares the bees' frequency when they are buzzing on 
# heavy (metal added) vs. light flowers
# Note: all pores were glued shut
#install packages
ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects", "gamm4", "viridis")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Loading required package: gamm4
## Loading required package: mgcv
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
## This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.
## This is gamm4 0.2-5
## Loading required package: viridis
## Loading required package: viridisLite
##  ggplot2 reshape2     lme4   sjPlot multcomp     plyr  effects    gamm4 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE 
##  viridis 
##     TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))

# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"


# print system info
print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run  2017-12-14 18:37:52"
print(R.version)
##                _                           
## platform       x86_64-apple-darwin15.6.0   
## arch           x86_64                      
## os             darwin15.6.0                
## system         x86_64, darwin15.6.0        
## status                                     
## major          3                           
## minor          4.2                         
## year           2017                        
## month          09                          
## day            28                          
## svn rev        73368                       
## language       R                           
## version.string R version 3.4.2 (2017-09-28)
## nickname       Short Summer

Read in data and do some visualizations

new_df = read.csv(file.path(dataDir, "02_HeavyLight_cleaned.csv"))
new_df$hive = as.factor(new_df$hive)
new_df$treatment <- mapvalues(new_df$treatment, from = c("sham", "weighted"), to = c("Sham", "Weighted"))
head(new_df)
##   freq     amp                  datetime rewNum
## 1  250 0.17945  2016_09_28__09_59_13_112      1
## 2  300 0.56657  2016_09_28__09_59_14_520      2
## 3  400 0.68371  2016_09_28__09_59_15_030      3
## 4  410 0.80118  2016_09_28__09_59_15_707      4
## 5  370 0.52666  2016_09_28__09_59_16_277      5
## 6  360 0.31541  2016_09_28__09_59_17_096      6
##                                     accFile beeID hive   IT
## 1 2016_09_28__09_59_13_112_220_450_test.txt     7    3 3.98
## 2 2016_09_28__09_59_14_520_220_450_test.txt     7    3 3.98
## 3 2016_09_28__09_59_15_030_220_450_test.txt     7    3 3.98
## 4 2016_09_28__09_59_15_707_220_450_test.txt     7    3 3.98
## 5 2016_09_28__09_59_16_277_220_450_test.txt     7    3 3.98
## 6 2016_09_28__09_59_17_096_220_450_test.txt     7    3 3.98
##                                             accFileAndFolder
## 1 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 2 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 3 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 4 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 5 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
## 6 Bee7_28Sept2016_Hive3_S_W/2016_09_28__09_58_58_ampFreq.txt
##   dateFrozenOrMarked treatment  amp_acc
## 1          28-Sep-16      Sham 17.64503
## 2          28-Sep-16      Sham 55.70993
## 3          28-Sep-16      Sham 67.22812
## 4          28-Sep-16      Sham 78.77876
## 5          28-Sep-16      Sham 51.78564
## 6          28-Sep-16      Sham 31.01377
# calculate average number of each treatment
ll = list(tapply(new_df$treatment, new_df$beeID, table))
dd = data.frame(unlist(ll))
dd$beeTrt = row.names(dd)
dd$beeID = sapply(dd$beeTrt, FUN = function(x) strsplit(x, "\\.")[[1]][1])
dd$trt = sapply(dd$beeTrt, FUN = function(x) strsplit(x, "\\.")[[1]][2])

# shows average number of buzzes per treatment
tapply(dd$unlist.ll., INDEX = dd$trt, mean)
##     Sham Weighted 
## 29.69444 35.86111
tapply(dd$unlist.ll., INDEX = dd$trt, sd) # sd for number of buzzes
##      Sham  Weighted 
##  9.603034 15.224328
# show histograms of frequencies
hist(new_df$freq, breaks = seq(215, 455, by = 10))

# Calculate average for each bee and for each treatment
frqMeans <- as.data.frame(tapply(X = new_df$freq, INDEX = list(new_df$beeID, new_df$treatment), mean))
frqMeans$beeID <- row.names(frqMeans)
frqMeans
##            Sham Weighted  beeID
## 1      344.3333 347.0833      1
## 10     331.2000 382.5000     10
## 11     411.7391 387.6923     11
## 12     344.8276 339.2500     12
## 14     324.0625 281.6000     14
## 15     367.5000 376.3158     15
## 16     357.0588 386.0606     16
## 17     314.0741 341.6000     17
## 18     366.4706 345.0000     18
## 19     375.0000 379.2308     19
## 2      298.5294 322.9167      2
## 20     359.8148 313.3333     20
## 21     331.3158 310.7143     21
## 22     399.4286 388.2857     22
## 23     348.2759 351.2821     23
## 24     287.7778 273.3333     24
## 25     345.0000 352.3684     25
## 26     369.6000 392.7273     26
## 27     366.6667 350.0000     27
## 28     364.4444 374.1935     28
## 29     291.2500 319.3333     29
## 3      399.0323 388.9655      3
## 30     368.0769 369.0323     30
## 31     273.5294 320.4167     31
## 32     325.6000 246.1538     32
## 4      356.6667 371.3636      4
## 5      370.7692 353.7143      5
## 6      284.2857 326.5714      6
## 7      366.9697 396.3415      7
## 8      333.8462 342.0833      8
## 9      344.4000 384.4444      9
## blue   319.1429 333.8462   blue
## gold   374.5098 365.2174   gold
## orange 309.6875 319.1000 orange
## pink   385.3333 367.5000   pink
## white  319.7674 326.9444  white
# convert to long format for ggplot-ing
frqLong <- melt(frqMeans, id.vars = "beeID", measure.vars =  c("Sham", "Weighted"), 
     variable.name = "trt", value.name = "frq")
nrow(frqLong)
## [1] 72
ggplot(frqLong, aes(x = trt, y = frq)) + 
     geom_boxplot() + 
    geom_line(aes(group = beeID))+
     geom_point()

ggplot(frqLong, aes(x = trt, y = frq)) + 
     geom_boxplot() + 
     labs(y = "Buzz Frequency (Hz)", x = "Flower treatment")

#ggsave(file.path(figDir, "heavyLight_frq.pdf"), width = 5, height = 4)

# print some descriptive statistics
nrow(new_df)
## [1] 2360
unique(new_df$hive)
## [1] 3 4
## Levels: 3 4
# Calculate average amplitude for each bee and for each treatment
ampMeans <- as.data.frame(tapply(X = new_df$amp_acc, INDEX = list(new_df$beeID, new_df$treatment), mean))

ampMeans$beeID <- row.names(ampMeans)
ampMeans
##            Sham  Weighted  beeID
## 1      79.41488 49.082432      1
## 10     48.45538 30.837297     10
## 11     34.56381 22.888359     11
## 12     39.69695 36.219420     12
## 14     28.06241 22.151111     14
## 15     43.84567 39.297029     15
## 16     50.08044 36.654450     16
## 17     38.58280 34.215988     17
## 18     31.56744 23.471063     18
## 19     43.66984 24.653556     19
## 2      35.11669 34.075563      2
## 20     47.39991 35.394684     20
## 21     44.71363 34.498420     21
## 22     32.36932 24.702542     22
## 23     63.43800 46.047727     23
## 24     41.36935 22.975869     24
## 25     36.78398 20.052036     25
## 26     28.09219 22.162838     26
## 27     39.66976 30.902595     27
## 28     44.74322 39.945697     28
## 29     43.33800 22.381776     29
## 3      49.02683 56.246838      3
## 30     54.47406 29.037523     30
## 31     57.17575 30.517863     31
## 32     35.02128  5.009379     32
## 4      34.84784 37.962188      4
## 5      51.37301 32.967355      5
## 6      25.06068 36.736087      6
## 7      42.30789 47.528431      7
## 8      71.53457 28.270977      8
## 9      38.79524 29.125278      9
## blue   31.83149 19.595870   blue
## gold   34.72996 20.867171   gold
## orange 52.83226 46.007837 orange
## pink   30.62350 39.548159   pink
## white  22.55167 15.381023  white
ampLong <- melt(ampMeans, id.vars = "beeID", measure.vars =  c("Sham", "Weighted"), 
                variable.name = "trt", value.name = "amp")

ggplot(ampLong, aes(x = trt, y = amp)) + 
     geom_boxplot() + 
     geom_point()  + 
  labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment")

ggplot(ampLong, aes(x = trt, y = amp)) + 
     geom_boxplot() + 
     geom_point() +
  geom_line(aes(group = beeID))

ggplot(ampLong, aes(x = trt, y = amp)) + 
     geom_boxplot() + 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment")

#ggsave(file.path(figDir, 'heavyLightAmp.pdf'), width = 5, heigh = 4)

nrow(ampLong)
## [1] 72

Modeling frequency with GAMM

Use BIC to select model (decide what interactions and covariates to use)

Make sure that REML = FALSE when comparing BIC values

Start with the biggest model of interest, and then see what predictors can be removed

rr = new_df$treatment[1]
df1 = new_df[1, ]
for(ii in 1:(nrow(new_df) -1)){
  if(rr == new_df$treatment[ii + 1]) next
  else{
    df1  = rbind(df1, new_df[ii+1, ])
    rr = new_df$treatment[ii + 1]
  }
}

df1 <- df1[df1$rewNum != 1, ]
df1$trtSwitch = df1$rewNum
colnames(df1)
##  [1] "freq"               "amp"                "datetime"          
##  [4] "rewNum"             "accFile"            "beeID"             
##  [7] "hive"               "IT"                 "accFileAndFolder"  
## [10] "dateFrozenOrMarked" "treatment"          "amp_acc"           
## [13] "trtSwitch"
df2 <- df1[, c("beeID", "hive", "trtSwitch", "IT")]
df2
##       beeID hive trtSwitch   IT
## 34        7    3        34 3.98
## 111   white    3        37 3.70
## 178       1    4        25 4.07
## 232       8    4        25 3.84
## 284      19    3        27 4.01
## 353      29    3        31 3.78
## 402      16    3        34 4.11
## 478      26    3        26 3.72
## 536      32    3        26 3.72
## 573      15    4        25 3.61
## 635      24    4        25 3.55
## 687      10    3        26 4.09
## 745      30    4        27 4.19
## 814      21    3        39 3.52
## 897      28    4        28 3.57
## 963       6    3        36 4.55
## 1002 orange    3        33 4.81
## 1129     17    4        28 4.25
## 1205     25    3        27 3.73
## 1267     31    3        25 3.76
## 1336   blue    3        36 3.91
## 1391     23    3        30 4.18
## 1459     12    3        30 4.10
## 1530      3    3        32 4.27
## 1594     22    4        36 3.65
## 1673      4    3        45 4.27
## 1710      9    3        26 4.17
## 1763      5    3        27 3.76
## 1832      2    4        35 4.42
## 1930     11    3        27 3.92
## 1986     20    3        34 4.28
## 2065     14    4        26 4.23
## 2127   pink    4        31 4.29
## 2204     27    3        34 3.51
## 2253     18    3        29 3.12
## 2310   gold    3        24 4.13
sum(duplicated(df2$beeID))
## [1] 0
newdf2 <- merge(new_df, df2)
head(newdf2)
##   beeID hive   IT freq     amp                  datetime rewNum
## 1     1    4 4.07  310 0.65224  2016_09_27__10_52_32_696      1
## 2     1    4 4.07  330 0.47733  2016_09_27__10_52_41_293      2
## 3     1    4 4.07  240 0.55067  2016_09_27__10_52_42_247      3
## 4     1    4 4.07  290 0.28977  2016_09_27__10_52_42_974      4
## 5     1    4 4.07  360 0.47311  2016_09_27__10_52_51_809      5
## 6     1    4 4.07  360 0.45364  2016_09_27__10_53_00_124      6
##                                     accFile
## 1 2016_09_27__10_52_32_696_220_450_test.txt
## 2 2016_09_27__10_52_41_293_220_450_test.txt
## 3 2016_09_27__10_52_42_247_220_450_test.txt
## 4 2016_09_27__10_52_42_974_220_450_test.txt
## 5 2016_09_27__10_52_51_809_220_450_test.txt
## 6 2016_09_27__10_53_00_124_220_450_test.txt
##                                         accFileAndFolder
## 1 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 2 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 3 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 4 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 5 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 6 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
##   dateFrozenOrMarked treatment  amp_acc trtSwitch
## 1          27-Sep-16  Weighted 64.13373        25
## 2          27-Sep-16  Weighted 46.93510        25
## 3          27-Sep-16  Weighted 54.14651        25
## 4          27-Sep-16  Weighted 28.49263        25
## 5          27-Sep-16  Weighted 46.52016        25
## 6          27-Sep-16  Weighted 44.60570        25
newdf2$visitNum_centered <- newdf2$rewNum - newdf2$trtSwitch

# get treatment
trts <- (sapply(unique(newdf2$beeID), FUN = function(x){
  startt = newdf2[newdf2$beeID == x & newdf2$rewNum == 1, "treatment" ]
  return(as.character(startt))
}))

df4 <- data.frame(beeID = unique(newdf2$beeID), trts = trts)

newdf3 <- merge(newdf2, df4)

newdf3$trt2 <- as.factor(paste(newdf3$trts, "first"))

head(newdf3)
##   beeID hive   IT freq     amp                  datetime rewNum
## 1     1    4 4.07  310 0.65224  2016_09_27__10_52_32_696      1
## 2     1    4 4.07  330 0.47733  2016_09_27__10_52_41_293      2
## 3     1    4 4.07  240 0.55067  2016_09_27__10_52_42_247      3
## 4     1    4 4.07  290 0.28977  2016_09_27__10_52_42_974      4
## 5     1    4 4.07  360 0.47311  2016_09_27__10_52_51_809      5
## 6     1    4 4.07  360 0.45364  2016_09_27__10_53_00_124      6
##                                     accFile
## 1 2016_09_27__10_52_32_696_220_450_test.txt
## 2 2016_09_27__10_52_41_293_220_450_test.txt
## 3 2016_09_27__10_52_42_247_220_450_test.txt
## 4 2016_09_27__10_52_42_974_220_450_test.txt
## 5 2016_09_27__10_52_51_809_220_450_test.txt
## 6 2016_09_27__10_53_00_124_220_450_test.txt
##                                         accFileAndFolder
## 1 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 2 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 3 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 4 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 5 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
## 6 Bee1_27Sept_Hive4_W_S/2016_09_27__10_52_05_ampFreq.txt
##   dateFrozenOrMarked treatment  amp_acc trtSwitch visitNum_centered
## 1          27-Sep-16  Weighted 64.13373        25               -24
## 2          27-Sep-16  Weighted 46.93510        25               -23
## 3          27-Sep-16  Weighted 54.14651        25               -22
## 4          27-Sep-16  Weighted 28.49263        25               -21
## 5          27-Sep-16  Weighted 46.52016        25               -20
## 6          27-Sep-16  Weighted 44.60570        25               -19
##       trts           trt2
## 1 Weighted Weighted first
## 2 Weighted Weighted first
## 3 Weighted Weighted first
## 4 Weighted Weighted first
## 5 Weighted Weighted first
## 6 Weighted Weighted first
# GAMM
# start with gamm so I can show change by visit number
g01 = gamm4(freq ~ s(visitNum_centered, by = trt2) + IT + hive + treatment + trt2, random =  ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)

summary(g01$gam) # Summary for paper 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## freq ~ s(visitNum_centered, by = trt2) + IT + hive + treatment + 
##     trt2
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         373.350     64.000   5.834 6.17e-09 ***
## IT                   -5.075     15.811  -0.321    0.748    
## hive4                -4.208     11.433  -0.368    0.713    
## treatmentWeighted     1.919      3.431   0.559    0.576    
## trt2Weighted first  -10.607     10.771  -0.985    0.325    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                         edf Ref.df     F p-value
## s(visitNum_centered):trt2Sham first       1      1 1.470   0.225
## s(visitNum_centered):trt2Weighted first   1      1 1.015   0.314
## 
## R-sq.(adj) =  -0.000401   
## lmer.REML =  25097  Scale est. = 2354.9    n = 2360
summary(g01$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 25096.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3193 -0.4698  0.1892  0.6934  2.3814 
## 
## Random effects:
##  Groups   Name                                    Variance  Std.Dev. 
##  beeID    (Intercept)                             9.602e+02 3.099e+01
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.000e+00 0.000e+00
##  Xr       s(visitNum_centered):trt2Sham first     3.028e-11 5.503e-06
##  Residual                                         2.355e+03 4.853e+01
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 373.350     64.000   5.834
## XIT                                           -5.075     15.811  -0.321
## Xhive4                                        -4.208     11.433  -0.368
## XtreatmentWeighted                             1.919      3.431   0.559
## Xtrt2Weighted first                          -10.607     10.771  -0.985
## Xs(visitNum_centered):trt2Sham firstFx1        2.285      1.885   1.212
## Xs(visitNum_centered):trt2Weighted firstFx1    2.569      2.550   1.007
## 
## Correlation of Fixed Effects:
##             X(Int) XIT    Xhive4 XtrtmW Xtr2Wf X(N_):2Sf
## XIT         -0.992                                      
## Xhive4      -0.043 -0.011                               
## XtrtmntWght -0.041  0.011  0.003                        
## Xtrt2Wghtdf -0.237  0.165 -0.014  0.051                 
## X(N_):2SfF1  0.049 -0.028 -0.008 -0.707 -0.038          
## X(N_):2WfF1 -0.032  0.012 -0.002  0.658  0.063 -0.466
dev.off()
## null device 
##           1
g02 = gamm4(freq ~ s(visitNum_centered, by = trt2) +  hive + treatment + trt2, random =  ~(1|beeID), data = newdf3)
anova(g01$mer, g02$mer) #g02 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g02$mer: NULL
## g01$mer: NULL
##         Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## g02$mer 10 25152 25210 -12566    25132                         
## g01$mer 11 25154 25218 -12566    25132 0.1168      1     0.7325
summary(g02$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 25104.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3173 -0.4693  0.1882  0.6954  2.3808 
## 
## Random effects:
##  Groups   Name                                    Variance  Std.Dev. 
##  beeID    (Intercept)                             9.331e+02 3.055e+01
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.000e+00 0.000e+00
##  Xr       s(visitNum_centered):trt2Sham first     1.379e-10 1.174e-05
##  Residual                                         2.355e+03 4.853e+01
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 352.967      7.981   44.23
## Xhive4                                        -4.248     11.276   -0.38
## XtreatmentWeighted                             1.939      3.431    0.57
## Xtrt2Weighted first                          -10.030     10.478   -0.96
## Xs(visitNum_centered):trt2Sham firstFx1        2.262      1.884    1.20
## Xs(visitNum_centered):trt2Weighted firstFx1    2.583      2.550    1.01
## 
## Correlation of Fixed Effects:
##             X(Int) Xhive4 XtrtmW Xtr2Wf X(N_):2Sf
## Xhive4      -0.426                               
## XtrtmntWght -0.242  0.004                        
## Xtrt2Wghtdf -0.587 -0.012  0.051                 
## X(N_):2SfF1  0.171 -0.009 -0.707 -0.034          
## X(N_):2WfF1 -0.158 -0.002  0.658  0.063 -0.466
g03 =  gamm4(freq ~ s(visitNum_centered, by = trt2) +  hive + trt2, random =  ~(1|beeID), data = newdf3)
anova(g02$mer, g03$mer) #g03 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g03$mer: NULL
## g02$mer: NULL
##         Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## g03$mer  9 25151 25203 -12566    25133                        
## g02$mer 10 25152 25210 -12566    25132 0.328      1     0.5668
summary(g03$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 25109
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3127 -0.4673  0.1923  0.6924  2.3917 
## 
## Random effects:
##  Groups   Name                                    Variance Std.Dev.
##  beeID    (Intercept)                              935.9   30.59   
##  Xr.0     s(visitNum_centered):trt2Weighted first    0.0    0.00   
##  Xr       s(visitNum_centered):trt2Sham first        0.0    0.00   
##  Residual                                         2354.1   48.52   
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 354.059      7.754   45.66
## Xhive4                                        -4.271     11.292   -0.38
## Xtrt2Weighted first                          -10.330     10.479   -0.99
## Xs(visitNum_centered):trt2Sham firstFx1        3.015      1.332    2.26
## Xs(visitNum_centered):trt2Weighted firstFx1    1.635      1.919    0.85
## 
## Correlation of Fixed Effects:
##             X(Int) Xhive4 Xtr2Wf X(N_):2Sf
## Xhive4      -0.438                        
## Xtrt2Wghtdf -0.593 -0.012                 
## X(N_):2SfF1  0.000 -0.008  0.003          
## X(N_):2WfF1  0.002 -0.005  0.039  0.000
g04 = gamm4(freq ~ s(visitNum_centered, by = trt2) +trt2, random =  ~(1|beeID), data = newdf3)
anova(g03$mer, g04$mer) # g04 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g04$mer: NULL
## g03$mer: NULL
##         Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## g04$mer  8 25149 25195 -12566    25133                         
## g03$mer  9 25151 25203 -12566    25133 0.1558      1      0.693
summary(g04$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 25115.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3139 -0.4690  0.1912  0.6937  2.3917 
## 
## Random effects:
##  Groups   Name                                    Variance  Std.Dev. 
##  beeID    (Intercept)                             9.112e+02 3.019e+01
##  Xr.0     s(visitNum_centered):trt2Weighted first 1.513e-09 3.890e-05
##  Xr       s(visitNum_centered):trt2Sham first     3.322e-12 1.823e-06
##  Residual                                         2.354e+03 4.852e+01
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 352.774      6.883   51.25
## Xtrt2Weighted first                          -10.376     10.345   -1.00
## Xs(visitNum_centered):trt2Sham firstFx1        3.008      1.332    2.26
## Xs(visitNum_centered):trt2Weighted firstFx1    1.631      1.919    0.85
## 
## Correlation of Fixed Effects:
##             X(Int) Xtr2Wf X(N_):2Sf
## Xtrt2Wghtdf -0.665                 
## X(N_):2SfF1 -0.004  0.003          
## X(N_):2WfF1  0.000  0.039  0.000
g05 = gamm4(freq ~ s(visitNum_centered) +trt2, random =  ~(1|beeID), data = newdf3)
anova(g04$mer, g05$mer) # g05 slightly better (according to BIC)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g05$mer: NULL
## g04$mer: NULL
##         Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## g05$mer  6 25145 25180 -12567    25133                         
## g04$mer  8 25149 25195 -12566    25133 0.3437      2     0.8421
summary(g05$gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## freq ~ s(visitNum_centered) + trt2
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         352.783      6.875  51.312   <2e-16 ***
## trt2Weighted first  -10.188     10.329  -0.986    0.324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df     F p-value  
## s(visitNum_centered)   1      1 5.478  0.0193 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  -0.00388   
## lmer.REML =  25120  Scale est. = 2353.5    n = 2360
g06 = gamm4(freq ~ s(visitNum_centered), random =  ~(1|beeID), data = newdf3)
anova(g05$mer, g06$mer) # g05 slightly better (according to BIC)
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g06$mer: NULL
## g05$mer: NULL
##         Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## g06$mer  5 25144 25173 -12567    25134                         
## g05$mer  6 25145 25180 -12567    25133 1.0153      1     0.3136
summary(g06$mer) # g06 Better
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 25127.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3036 -0.4740  0.1917  0.6901  2.3998 
## 
## Random effects:
##  Groups   Name                 Variance  Std.Dev. 
##  beeID    (Intercept)          9.078e+02 3.013e+01
##  Xr       s(visitNum_centered) 9.414e-11 9.703e-06
##  Residual                      2.354e+03 4.851e+01
## Number of obs: 2360, groups:  beeID, 36; Xr, 8
## 
## Fixed effects:
##                          Estimate Std. Error t value
## X(Intercept)              348.268      5.127   67.93
## Xs(visitNum_centered)Fx1    2.587      1.093    2.37
## 
## Correlation of Fixed Effects:
##             X(Int)
## Xs(vstN_)F1 0.018
g07 = gamm4(freq ~ 1, random =  ~(1|beeID), data = newdf3)
anova(g06$mer, g07$mer) # g07 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g07$mer: NULL
## g06$mer: NULL
##         Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## g07$mer  3 25146 25163 -12570    25140                          
## g06$mer  5 25144 25173 -12567    25134 5.584      2     0.0613 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# looks like GAMM is not necessary -- no predictors significantly affect frequency

Modeling frequency with lmer

Do a similar thing to GAMM above (but make it more interpretable)

# fit lmer that is similar to gamm above
# interaction that may be important, based on domain knowledge
m0 = lmer(freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive  + I(scale(visitNum_centered^2)) + treatment + (1|beeID), data = newdf3, REML = FALSE)

summary(m0)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive + I(scale(visitNum_centered^2)) +  
##     treatment + (1 | beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25152.3  25210.0 -12566.1  25132.3     2350 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3176 -0.4731  0.1881  0.6953  2.3769 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  848.9   29.14   
##  Residual             2351.9   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                                                Estimate Std. Error t value
## (Intercept)                                    373.6133    60.3485   6.191
## I(scale(visitNum_centered))                      2.0637     2.4466   0.844
## trt2Weighted first                             -10.5531    10.1598  -1.039
## IT                                              -5.1688    14.9133  -0.347
## hive4                                           -4.1822    10.7784  -0.388
## I(scale(visitNum_centered^2))                    0.1987     1.5795   0.126
## treatmentWeighted                                2.1164     3.6713   0.576
## I(scale(visitNum_centered)):trt2Weighted first   0.5702     4.2740   0.133
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_) trt2Wf IT     hive4  I((N_^ trtmnW
## I(scl(vN_))  0.023                                          
## trt2Wghtdfr -0.237 -0.049                                   
## IT          -0.992  0.002  0.164                            
## hive4       -0.042 -0.018 -0.013 -0.012                     
## I(sc(N_^2))  0.026 -0.638  0.029 -0.038  0.018              
## trtmntWghtd -0.031 -0.737  0.061 -0.003  0.010  0.357       
## I((N_)):2Wf -0.031 -0.842  0.070  0.004  0.011  0.453  0.820
# no interaction
m1 = update(m0, .~. - IT)
BIC(m0, m1) # stay with m1 -- no interaction
##    df      BIC
## m0 10 25209.95
## m1  9 25202.31
anova(m0, m1) # agrees with BIC
## Data: newdf3
## Models:
## m1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) + 
## m1:     treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
## m0: freq ~ I(scale(visitNum_centered)) * trt2 + IT + hive + I(scale(visitNum_centered^2)) + 
## m0:     treatment + (1 | beeID)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m1  9 25150 25202 -12566    25132                         
## m0 10 25152 25210 -12566    25132 0.1199      1     0.7291
summary(m1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) +  
##     treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25150.4  25202.3 -12566.2  25132.4     2351 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3166 -0.4731  0.1875  0.6966  2.3778 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  852.1   29.19   
##  Residual             2351.8   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                                                Estimate Std. Error t value
## (Intercept)                                    352.8677     7.6941   45.86
## I(scale(visitNum_centered))                      2.0655     2.4466    0.84
## trt2Weighted first                              -9.9744    10.0391   -0.99
## hive4                                           -4.2258    10.7968   -0.39
## I(scale(visitNum_centered^2))                    0.1783     1.5784    0.11
## treatmentWeighted                                2.1121     3.6713    0.58
## I(scale(visitNum_centered)):trt2Weighted first   0.5745     4.2740    0.13
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_) trt2Wf hive4  I((N_^ trtmnW
## I(scl(vN_))  0.196                                   
## trt2Wghtdfr -0.586 -0.050                            
## hive4       -0.424 -0.018 -0.011                     
## I(sc(N_^2)) -0.093 -0.639  0.036  0.018              
## trtmntWghtd -0.268 -0.737  0.062  0.010  0.358       
## I((N_)):2Wf -0.218 -0.842  0.071  0.011  0.454  0.820
m1.1 <- update(m1, .~. - I(scale(visitNum_centered^2)))
anova(m1.1, m1) # keep 1.1
## Data: newdf3
## Models:
## m1.1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment + 
## m1.1:     (1 | beeID) + I(scale(visitNum_centered)):trt2
## m1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + I(scale(visitNum_centered^2)) + 
## m1:     treatment + (1 | beeID) + I(scale(visitNum_centered)):trt2
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m1.1  8 25148 25195 -12566    25132                         
## m1    9 25150 25202 -12566    25132 0.0128      1     0.9101
summary(m1.1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +  
##     (1 | beeID) + I(scale(visitNum_centered)):trt2
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25148.4  25194.6 -12566.2  25132.4     2352 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3159 -0.4706  0.1876  0.6974  2.3783 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  851.6   29.18   
##  Residual             2351.9   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                                                Estimate Std. Error t value
## (Intercept)                                    352.9488     7.6583   46.09
## I(scale(visitNum_centered))                      2.2418     1.8826    1.19
## trt2Weighted first                             -10.0149    10.0298   -1.00
## hive4                                           -4.2477    10.7920   -0.39
## treatmentWeighted                                1.9640     3.4286    0.57
## I(scale(visitNum_centered)):trt2Weighted first   0.3557     3.8085    0.09
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_) trt2Wf hive4  trtmnW
## I(scl(vN_))  0.178                            
## trt2Wghtdfr -0.586 -0.035                     
## hive4       -0.425 -0.009 -0.012              
## trtmntWghtd -0.252 -0.707  0.053  0.004       
## I((N_)):2Wf -0.198 -0.806  0.061  0.003  0.790
m2 <- update(m1.1, .~. - I(scale(visitNum_centered)):trt2)
anova(m1.1, m2) # keep m2 (based on BIC)
## Data: newdf3
## Models:
## m2: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment + 
## m2:     (1 | beeID)
## m1.1: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment + 
## m1.1:     (1 | beeID) + I(scale(visitNum_centered)):trt2
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2    7 25146 25187 -12566    25132                         
## m1.1  8 25148 25195 -12566    25132 0.0087      1     0.9256
summary(m2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment +  
##     (1 | beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25146.4  25186.8 -12566.2  25132.4     2353 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3167 -0.4685  0.1872  0.6957  2.3784 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  852.1   29.19   
##  Residual             2351.9   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                  353.090      7.509   47.02
## I(scale(visitNum_centered))    2.384      1.115    2.14
## trt2Weighted first           -10.072     10.014   -1.01
## hive4                         -4.251     10.795   -0.39
## treatmentWeighted              1.711      2.102    0.81
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_) trt2Wf hive4 
## I(scl(vN_))  0.032                     
## trt2Wghtdfr -0.587  0.024              
## hive4       -0.433 -0.011 -0.012       
## trtmntWghtd -0.159 -0.194  0.007  0.002
m3 <- update(m2, .~. - hive)
anova(m2, m3) # keep m3
## Data: newdf3
## Models:
## m3: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 | 
## m3:     beeID)
## m2: freq ~ I(scale(visitNum_centered)) + trt2 + hive + treatment + 
## m2:     (1 | beeID)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  6 25145 25179 -12566    25133                         
## m2  7 25146 25187 -12566    25132 0.1547      1     0.6941
summary(m3)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 |  
##     beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25144.6  25179.2 -12566.3  25132.6     2354 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3192 -0.4693  0.1879  0.6954  2.3798 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  855.9   29.26   
##  Residual             2351.9   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                  351.811      6.784   51.86
## I(scale(visitNum_centered))    2.379      1.115    2.13
## trt2Weighted first           -10.121     10.035   -1.01
## treatmentWeighted              1.713      2.102    0.81
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_) trt2Wf
## I(scl(vN_))  0.031              
## trt2Wghtdfr -0.657  0.024       
## trtmntWghtd -0.175 -0.194  0.007
m3.1 <- update(m3, .~. - treatment)
anova(m3, m3.1) # keep m3.1
## Data: newdf3
## Models:
## m3.1: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
## m3: freq ~ I(scale(visitNum_centered)) + trt2 + treatment + (1 | 
## m3:     beeID)
##      Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3.1  5 25143 25172 -12567    25133                        
## m3    6 25145 25179 -12566    25133 0.664      1     0.4151
summary(m3.1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25143.2  25172.1 -12566.6  25133.2     2355 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2984 -0.4723  0.1895  0.6903  2.4024 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  856.2   29.26   
##  Residual             2352.5   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                  352.780      6.680   52.81
## I(scale(visitNum_centered))    2.556      1.094    2.34
## trt2Weighted first           -10.181     10.036   -1.01
## 
## Correlation of Fixed Effects:
##             (Intr) I((N_)
## I(scl(vN_)) -0.003       
## trt2Wghtdfr -0.666  0.025
m4 <- update(m3.1, .~. - trt2)
summary(m4)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##  25142.3  25165.3 -12567.1  25134.3     2356 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3033 -0.4735  0.1920  0.6906  2.3989 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  881.3   29.69   
##  Residual             2352.5   48.50   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                             Estimate Std. Error t value
## (Intercept)                  348.269      5.054   68.91
## I(scale(visitNum_centered))    2.586      1.093    2.36
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(scl(vN_)) 0.018
anova(m3.1, m4) #keep m4
## Data: newdf3
## Models:
## m4: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
## m3.1: freq ~ I(scale(visitNum_centered)) + trt2 + (1 | beeID)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m4    4 25142 25165 -12567    25134                         
## m3.1  5 25143 25172 -12567    25133 1.0153      1     0.3136
m5 <- update(m4, .~. -I(scale(visitNum_centered)))
anova(m4, m5) # keep m5 -- baesd on BIC (agrees with GAMM analysis above)
## Data: newdf3
## Models:
## m5: freq ~ (1 | beeID)
## m4: freq ~ I(scale(visitNum_centered)) + (1 | beeID)
##    Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## m5  3 25146 25163 -12570    25140                          
## m4  4 25142 25165 -12567    25134 5.584      1    0.01813 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# final freq mod for paper
finmod <- update(m5, .~., REML = TRUE)
summary(finmod) # final model for paper -- same as gamm
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID)
##    Data: newdf3
## 
## REML criterion at convergence: 25134.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3401 -0.4580  0.1806  0.6919  2.3616 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  901.6   30.03   
##  Residual             2358.4   48.56   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  348.056      5.109   68.12
summary(finmod) # summary for paper -- no predictors predict frequency (when accounting for beeID)
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID)
##    Data: newdf3
## 
## REML criterion at convergence: 25134.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3401 -0.4580  0.1806  0.6919  2.3616 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  901.6   30.03   
##  Residual             2358.4   48.56   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  348.056      5.109   68.12
plot(finmod)

qqnorm(ranef(finmod)$beeID[[1]])
qqline(ranef(finmod)$beeID[[1]])

sjp.lmer(finmod, type = "re", sort = TRUE) # plot random effects to find any outliers
## Plotting random effects...

# sjp.lmer(finmod, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects

Bootstrap CI’s for figure for paper

# set number of bootstrap replicates for models
nbootSims = 1000

table(new_df$hive) # more trials from hive 3
## 
##    3    4 
## 1644  716
# using hive 3, since it's the one with the most data
# however, hive doesn't affect model anyway

# calculate an average IT for prediction
ITmean = mean(tapply(new_df$IT, INDEX = new_df$beeID, FUN = function(x) x[1] ))

pframe <- expand.grid(trt2 = levels(droplevels(newdf3$trt2)),
                     IT = ITmean,
                     visitNum_centered = -50:50, 
                     hive = factor(3, levels = levels(newdf3$hive)),  
                     beeID = 99999, 
                     treatment = levels(newdf3$treatment))
pframe$freq <- 0
pframe <- pframe[(pframe$trt2 == "Sham first" &
                    pframe$treatment == "Sham" & 
                    pframe$visitNum_centered < 0) |
                    (pframe$trt2 == "Sham first" &
                    pframe$treatment == "Weighted" & 
                    pframe$visitNum_centered > 0) |
                    (pframe$trt2 == "Weighted first" &
                    pframe$treatment == "Weighted" & 
                    pframe$visitNum_centered < 0) |
                    (pframe$trt2 == "Weighted first" &
                    pframe$treatment == "Sham" & 
                    pframe$visitNum_centered > 0) 
                    , ]

# make a model to show the effect of treatment in a plot
m1 <- update(finmod, .~. + treatment)
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
## Warning in optwrap(object@optinfo$optimizer, ff, x0, lower = lower, control
## = control$optCtrl, : convergence code 3 from bobyqa: bobyqa -- a trust
## region step failed to reduce q
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe
##               trt2       IT visitNum_centered hive beeID treatment freq
## 1       Sham first 3.965833               -50    3 99999      Sham    0
## 3       Sham first 3.965833               -49    3 99999      Sham    0
## 5       Sham first 3.965833               -48    3 99999      Sham    0
## 7       Sham first 3.965833               -47    3 99999      Sham    0
## 9       Sham first 3.965833               -46    3 99999      Sham    0
## 11      Sham first 3.965833               -45    3 99999      Sham    0
## 13      Sham first 3.965833               -44    3 99999      Sham    0
## 15      Sham first 3.965833               -43    3 99999      Sham    0
## 17      Sham first 3.965833               -42    3 99999      Sham    0
## 19      Sham first 3.965833               -41    3 99999      Sham    0
## 21      Sham first 3.965833               -40    3 99999      Sham    0
## 23      Sham first 3.965833               -39    3 99999      Sham    0
## 25      Sham first 3.965833               -38    3 99999      Sham    0
## 27      Sham first 3.965833               -37    3 99999      Sham    0
## 29      Sham first 3.965833               -36    3 99999      Sham    0
## 31      Sham first 3.965833               -35    3 99999      Sham    0
## 33      Sham first 3.965833               -34    3 99999      Sham    0
## 35      Sham first 3.965833               -33    3 99999      Sham    0
## 37      Sham first 3.965833               -32    3 99999      Sham    0
## 39      Sham first 3.965833               -31    3 99999      Sham    0
## 41      Sham first 3.965833               -30    3 99999      Sham    0
## 43      Sham first 3.965833               -29    3 99999      Sham    0
## 45      Sham first 3.965833               -28    3 99999      Sham    0
## 47      Sham first 3.965833               -27    3 99999      Sham    0
## 49      Sham first 3.965833               -26    3 99999      Sham    0
## 51      Sham first 3.965833               -25    3 99999      Sham    0
## 53      Sham first 3.965833               -24    3 99999      Sham    0
## 55      Sham first 3.965833               -23    3 99999      Sham    0
## 57      Sham first 3.965833               -22    3 99999      Sham    0
## 59      Sham first 3.965833               -21    3 99999      Sham    0
## 61      Sham first 3.965833               -20    3 99999      Sham    0
## 63      Sham first 3.965833               -19    3 99999      Sham    0
## 65      Sham first 3.965833               -18    3 99999      Sham    0
## 67      Sham first 3.965833               -17    3 99999      Sham    0
## 69      Sham first 3.965833               -16    3 99999      Sham    0
## 71      Sham first 3.965833               -15    3 99999      Sham    0
## 73      Sham first 3.965833               -14    3 99999      Sham    0
## 75      Sham first 3.965833               -13    3 99999      Sham    0
## 77      Sham first 3.965833               -12    3 99999      Sham    0
## 79      Sham first 3.965833               -11    3 99999      Sham    0
## 81      Sham first 3.965833               -10    3 99999      Sham    0
## 83      Sham first 3.965833                -9    3 99999      Sham    0
## 85      Sham first 3.965833                -8    3 99999      Sham    0
## 87      Sham first 3.965833                -7    3 99999      Sham    0
## 89      Sham first 3.965833                -6    3 99999      Sham    0
## 91      Sham first 3.965833                -5    3 99999      Sham    0
## 93      Sham first 3.965833                -4    3 99999      Sham    0
## 95      Sham first 3.965833                -3    3 99999      Sham    0
## 97      Sham first 3.965833                -2    3 99999      Sham    0
## 99      Sham first 3.965833                -1    3 99999      Sham    0
## 104 Weighted first 3.965833                 1    3 99999      Sham    0
## 106 Weighted first 3.965833                 2    3 99999      Sham    0
## 108 Weighted first 3.965833                 3    3 99999      Sham    0
## 110 Weighted first 3.965833                 4    3 99999      Sham    0
## 112 Weighted first 3.965833                 5    3 99999      Sham    0
## 114 Weighted first 3.965833                 6    3 99999      Sham    0
## 116 Weighted first 3.965833                 7    3 99999      Sham    0
## 118 Weighted first 3.965833                 8    3 99999      Sham    0
## 120 Weighted first 3.965833                 9    3 99999      Sham    0
## 122 Weighted first 3.965833                10    3 99999      Sham    0
## 124 Weighted first 3.965833                11    3 99999      Sham    0
## 126 Weighted first 3.965833                12    3 99999      Sham    0
## 128 Weighted first 3.965833                13    3 99999      Sham    0
## 130 Weighted first 3.965833                14    3 99999      Sham    0
## 132 Weighted first 3.965833                15    3 99999      Sham    0
## 134 Weighted first 3.965833                16    3 99999      Sham    0
## 136 Weighted first 3.965833                17    3 99999      Sham    0
## 138 Weighted first 3.965833                18    3 99999      Sham    0
## 140 Weighted first 3.965833                19    3 99999      Sham    0
## 142 Weighted first 3.965833                20    3 99999      Sham    0
## 144 Weighted first 3.965833                21    3 99999      Sham    0
## 146 Weighted first 3.965833                22    3 99999      Sham    0
## 148 Weighted first 3.965833                23    3 99999      Sham    0
## 150 Weighted first 3.965833                24    3 99999      Sham    0
## 152 Weighted first 3.965833                25    3 99999      Sham    0
## 154 Weighted first 3.965833                26    3 99999      Sham    0
## 156 Weighted first 3.965833                27    3 99999      Sham    0
## 158 Weighted first 3.965833                28    3 99999      Sham    0
## 160 Weighted first 3.965833                29    3 99999      Sham    0
## 162 Weighted first 3.965833                30    3 99999      Sham    0
## 164 Weighted first 3.965833                31    3 99999      Sham    0
## 166 Weighted first 3.965833                32    3 99999      Sham    0
## 168 Weighted first 3.965833                33    3 99999      Sham    0
## 170 Weighted first 3.965833                34    3 99999      Sham    0
## 172 Weighted first 3.965833                35    3 99999      Sham    0
## 174 Weighted first 3.965833                36    3 99999      Sham    0
## 176 Weighted first 3.965833                37    3 99999      Sham    0
## 178 Weighted first 3.965833                38    3 99999      Sham    0
## 180 Weighted first 3.965833                39    3 99999      Sham    0
## 182 Weighted first 3.965833                40    3 99999      Sham    0
## 184 Weighted first 3.965833                41    3 99999      Sham    0
## 186 Weighted first 3.965833                42    3 99999      Sham    0
## 188 Weighted first 3.965833                43    3 99999      Sham    0
## 190 Weighted first 3.965833                44    3 99999      Sham    0
## 192 Weighted first 3.965833                45    3 99999      Sham    0
## 194 Weighted first 3.965833                46    3 99999      Sham    0
## 196 Weighted first 3.965833                47    3 99999      Sham    0
## 198 Weighted first 3.965833                48    3 99999      Sham    0
## 200 Weighted first 3.965833                49    3 99999      Sham    0
## 202 Weighted first 3.965833                50    3 99999      Sham    0
## 204 Weighted first 3.965833               -50    3 99999  Weighted    0
## 206 Weighted first 3.965833               -49    3 99999  Weighted    0
## 208 Weighted first 3.965833               -48    3 99999  Weighted    0
## 210 Weighted first 3.965833               -47    3 99999  Weighted    0
## 212 Weighted first 3.965833               -46    3 99999  Weighted    0
## 214 Weighted first 3.965833               -45    3 99999  Weighted    0
## 216 Weighted first 3.965833               -44    3 99999  Weighted    0
## 218 Weighted first 3.965833               -43    3 99999  Weighted    0
## 220 Weighted first 3.965833               -42    3 99999  Weighted    0
## 222 Weighted first 3.965833               -41    3 99999  Weighted    0
## 224 Weighted first 3.965833               -40    3 99999  Weighted    0
## 226 Weighted first 3.965833               -39    3 99999  Weighted    0
## 228 Weighted first 3.965833               -38    3 99999  Weighted    0
## 230 Weighted first 3.965833               -37    3 99999  Weighted    0
## 232 Weighted first 3.965833               -36    3 99999  Weighted    0
## 234 Weighted first 3.965833               -35    3 99999  Weighted    0
## 236 Weighted first 3.965833               -34    3 99999  Weighted    0
## 238 Weighted first 3.965833               -33    3 99999  Weighted    0
## 240 Weighted first 3.965833               -32    3 99999  Weighted    0
## 242 Weighted first 3.965833               -31    3 99999  Weighted    0
## 244 Weighted first 3.965833               -30    3 99999  Weighted    0
## 246 Weighted first 3.965833               -29    3 99999  Weighted    0
## 248 Weighted first 3.965833               -28    3 99999  Weighted    0
## 250 Weighted first 3.965833               -27    3 99999  Weighted    0
## 252 Weighted first 3.965833               -26    3 99999  Weighted    0
## 254 Weighted first 3.965833               -25    3 99999  Weighted    0
## 256 Weighted first 3.965833               -24    3 99999  Weighted    0
## 258 Weighted first 3.965833               -23    3 99999  Weighted    0
## 260 Weighted first 3.965833               -22    3 99999  Weighted    0
## 262 Weighted first 3.965833               -21    3 99999  Weighted    0
## 264 Weighted first 3.965833               -20    3 99999  Weighted    0
## 266 Weighted first 3.965833               -19    3 99999  Weighted    0
## 268 Weighted first 3.965833               -18    3 99999  Weighted    0
## 270 Weighted first 3.965833               -17    3 99999  Weighted    0
## 272 Weighted first 3.965833               -16    3 99999  Weighted    0
## 274 Weighted first 3.965833               -15    3 99999  Weighted    0
## 276 Weighted first 3.965833               -14    3 99999  Weighted    0
## 278 Weighted first 3.965833               -13    3 99999  Weighted    0
## 280 Weighted first 3.965833               -12    3 99999  Weighted    0
## 282 Weighted first 3.965833               -11    3 99999  Weighted    0
## 284 Weighted first 3.965833               -10    3 99999  Weighted    0
## 286 Weighted first 3.965833                -9    3 99999  Weighted    0
## 288 Weighted first 3.965833                -8    3 99999  Weighted    0
## 290 Weighted first 3.965833                -7    3 99999  Weighted    0
## 292 Weighted first 3.965833                -6    3 99999  Weighted    0
## 294 Weighted first 3.965833                -5    3 99999  Weighted    0
## 296 Weighted first 3.965833                -4    3 99999  Weighted    0
## 298 Weighted first 3.965833                -3    3 99999  Weighted    0
## 300 Weighted first 3.965833                -2    3 99999  Weighted    0
## 302 Weighted first 3.965833                -1    3 99999  Weighted    0
## 305     Sham first 3.965833                 1    3 99999  Weighted    0
## 307     Sham first 3.965833                 2    3 99999  Weighted    0
## 309     Sham first 3.965833                 3    3 99999  Weighted    0
## 311     Sham first 3.965833                 4    3 99999  Weighted    0
## 313     Sham first 3.965833                 5    3 99999  Weighted    0
## 315     Sham first 3.965833                 6    3 99999  Weighted    0
## 317     Sham first 3.965833                 7    3 99999  Weighted    0
## 319     Sham first 3.965833                 8    3 99999  Weighted    0
## 321     Sham first 3.965833                 9    3 99999  Weighted    0
## 323     Sham first 3.965833                10    3 99999  Weighted    0
## 325     Sham first 3.965833                11    3 99999  Weighted    0
## 327     Sham first 3.965833                12    3 99999  Weighted    0
## 329     Sham first 3.965833                13    3 99999  Weighted    0
## 331     Sham first 3.965833                14    3 99999  Weighted    0
## 333     Sham first 3.965833                15    3 99999  Weighted    0
## 335     Sham first 3.965833                16    3 99999  Weighted    0
## 337     Sham first 3.965833                17    3 99999  Weighted    0
## 339     Sham first 3.965833                18    3 99999  Weighted    0
## 341     Sham first 3.965833                19    3 99999  Weighted    0
## 343     Sham first 3.965833                20    3 99999  Weighted    0
## 345     Sham first 3.965833                21    3 99999  Weighted    0
## 347     Sham first 3.965833                22    3 99999  Weighted    0
## 349     Sham first 3.965833                23    3 99999  Weighted    0
## 351     Sham first 3.965833                24    3 99999  Weighted    0
## 353     Sham first 3.965833                25    3 99999  Weighted    0
## 355     Sham first 3.965833                26    3 99999  Weighted    0
## 357     Sham first 3.965833                27    3 99999  Weighted    0
## 359     Sham first 3.965833                28    3 99999  Weighted    0
## 361     Sham first 3.965833                29    3 99999  Weighted    0
## 363     Sham first 3.965833                30    3 99999  Weighted    0
## 365     Sham first 3.965833                31    3 99999  Weighted    0
## 367     Sham first 3.965833                32    3 99999  Weighted    0
## 369     Sham first 3.965833                33    3 99999  Weighted    0
## 371     Sham first 3.965833                34    3 99999  Weighted    0
## 373     Sham first 3.965833                35    3 99999  Weighted    0
## 375     Sham first 3.965833                36    3 99999  Weighted    0
## 377     Sham first 3.965833                37    3 99999  Weighted    0
## 379     Sham first 3.965833                38    3 99999  Weighted    0
## 381     Sham first 3.965833                39    3 99999  Weighted    0
## 383     Sham first 3.965833                40    3 99999  Weighted    0
## 385     Sham first 3.965833                41    3 99999  Weighted    0
## 387     Sham first 3.965833                42    3 99999  Weighted    0
## 389     Sham first 3.965833                43    3 99999  Weighted    0
## 391     Sham first 3.965833                44    3 99999  Weighted    0
## 393     Sham first 3.965833                45    3 99999  Weighted    0
## 395     Sham first 3.965833                46    3 99999  Weighted    0
## 397     Sham first 3.965833                47    3 99999  Weighted    0
## 399     Sham first 3.965833                48    3 99999  Weighted    0
## 401     Sham first 3.965833                49    3 99999  Weighted    0
## 403     Sham first 3.965833                50    3 99999  Weighted    0
##          blo      bhi predMean
## 1   335.8801 357.9694 346.6421
## 3   335.8801 357.9694 346.6421
## 5   335.8801 357.9694 346.6421
## 7   335.8801 357.9694 346.6421
## 9   335.8801 357.9694 346.6421
## 11  335.8801 357.9694 346.6421
## 13  335.8801 357.9694 346.6421
## 15  335.8801 357.9694 346.6421
## 17  335.8801 357.9694 346.6421
## 19  335.8801 357.9694 346.6421
## 21  335.8801 357.9694 346.6421
## 23  335.8801 357.9694 346.6421
## 25  335.8801 357.9694 346.6421
## 27  335.8801 357.9694 346.6421
## 29  335.8801 357.9694 346.6421
## 31  335.8801 357.9694 346.6421
## 33  335.8801 357.9694 346.6421
## 35  335.8801 357.9694 346.6421
## 37  335.8801 357.9694 346.6421
## 39  335.8801 357.9694 346.6421
## 41  335.8801 357.9694 346.6421
## 43  335.8801 357.9694 346.6421
## 45  335.8801 357.9694 346.6421
## 47  335.8801 357.9694 346.6421
## 49  335.8801 357.9694 346.6421
## 51  335.8801 357.9694 346.6421
## 53  335.8801 357.9694 346.6421
## 55  335.8801 357.9694 346.6421
## 57  335.8801 357.9694 346.6421
## 59  335.8801 357.9694 346.6421
## 61  335.8801 357.9694 346.6421
## 63  335.8801 357.9694 346.6421
## 65  335.8801 357.9694 346.6421
## 67  335.8801 357.9694 346.6421
## 69  335.8801 357.9694 346.6421
## 71  335.8801 357.9694 346.6421
## 73  335.8801 357.9694 346.6421
## 75  335.8801 357.9694 346.6421
## 77  335.8801 357.9694 346.6421
## 79  335.8801 357.9694 346.6421
## 81  335.8801 357.9694 346.6421
## 83  335.8801 357.9694 346.6421
## 85  335.8801 357.9694 346.6421
## 87  335.8801 357.9694 346.6421
## 89  335.8801 357.9694 346.6421
## 91  335.8801 357.9694 346.6421
## 93  335.8801 357.9694 346.6421
## 95  335.8801 357.9694 346.6421
## 97  335.8801 357.9694 346.6421
## 99  335.8801 357.9694 346.6421
## 104 335.8801 357.9694 346.6421
## 106 335.8801 357.9694 346.6421
## 108 335.8801 357.9694 346.6421
## 110 335.8801 357.9694 346.6421
## 112 335.8801 357.9694 346.6421
## 114 335.8801 357.9694 346.6421
## 116 335.8801 357.9694 346.6421
## 118 335.8801 357.9694 346.6421
## 120 335.8801 357.9694 346.6421
## 122 335.8801 357.9694 346.6421
## 124 335.8801 357.9694 346.6421
## 126 335.8801 357.9694 346.6421
## 128 335.8801 357.9694 346.6421
## 130 335.8801 357.9694 346.6421
## 132 335.8801 357.9694 346.6421
## 134 335.8801 357.9694 346.6421
## 136 335.8801 357.9694 346.6421
## 138 335.8801 357.9694 346.6421
## 140 335.8801 357.9694 346.6421
## 142 335.8801 357.9694 346.6421
## 144 335.8801 357.9694 346.6421
## 146 335.8801 357.9694 346.6421
## 148 335.8801 357.9694 346.6421
## 150 335.8801 357.9694 346.6421
## 152 335.8801 357.9694 346.6421
## 154 335.8801 357.9694 346.6421
## 156 335.8801 357.9694 346.6421
## 158 335.8801 357.9694 346.6421
## 160 335.8801 357.9694 346.6421
## 162 335.8801 357.9694 346.6421
## 164 335.8801 357.9694 346.6421
## 166 335.8801 357.9694 346.6421
## 168 335.8801 357.9694 346.6421
## 170 335.8801 357.9694 346.6421
## 172 335.8801 357.9694 346.6421
## 174 335.8801 357.9694 346.6421
## 176 335.8801 357.9694 346.6421
## 178 335.8801 357.9694 346.6421
## 180 335.8801 357.9694 346.6421
## 182 335.8801 357.9694 346.6421
## 184 335.8801 357.9694 346.6421
## 186 335.8801 357.9694 346.6421
## 188 335.8801 357.9694 346.6421
## 190 335.8801 357.9694 346.6421
## 192 335.8801 357.9694 346.6421
## 194 335.8801 357.9694 346.6421
## 196 335.8801 357.9694 346.6421
## 198 335.8801 357.9694 346.6421
## 200 335.8801 357.9694 346.6421
## 202 335.8801 357.9694 346.6421
## 204 339.2423 359.9579 349.2520
## 206 339.2423 359.9579 349.2520
## 208 339.2423 359.9579 349.2520
## 210 339.2423 359.9579 349.2520
## 212 339.2423 359.9579 349.2520
## 214 339.2423 359.9579 349.2520
## 216 339.2423 359.9579 349.2520
## 218 339.2423 359.9579 349.2520
## 220 339.2423 359.9579 349.2520
## 222 339.2423 359.9579 349.2520
## 224 339.2423 359.9579 349.2520
## 226 339.2423 359.9579 349.2520
## 228 339.2423 359.9579 349.2520
## 230 339.2423 359.9579 349.2520
## 232 339.2423 359.9579 349.2520
## 234 339.2423 359.9579 349.2520
## 236 339.2423 359.9579 349.2520
## 238 339.2423 359.9579 349.2520
## 240 339.2423 359.9579 349.2520
## 242 339.2423 359.9579 349.2520
## 244 339.2423 359.9579 349.2520
## 246 339.2423 359.9579 349.2520
## 248 339.2423 359.9579 349.2520
## 250 339.2423 359.9579 349.2520
## 252 339.2423 359.9579 349.2520
## 254 339.2423 359.9579 349.2520
## 256 339.2423 359.9579 349.2520
## 258 339.2423 359.9579 349.2520
## 260 339.2423 359.9579 349.2520
## 262 339.2423 359.9579 349.2520
## 264 339.2423 359.9579 349.2520
## 266 339.2423 359.9579 349.2520
## 268 339.2423 359.9579 349.2520
## 270 339.2423 359.9579 349.2520
## 272 339.2423 359.9579 349.2520
## 274 339.2423 359.9579 349.2520
## 276 339.2423 359.9579 349.2520
## 278 339.2423 359.9579 349.2520
## 280 339.2423 359.9579 349.2520
## 282 339.2423 359.9579 349.2520
## 284 339.2423 359.9579 349.2520
## 286 339.2423 359.9579 349.2520
## 288 339.2423 359.9579 349.2520
## 290 339.2423 359.9579 349.2520
## 292 339.2423 359.9579 349.2520
## 294 339.2423 359.9579 349.2520
## 296 339.2423 359.9579 349.2520
## 298 339.2423 359.9579 349.2520
## 300 339.2423 359.9579 349.2520
## 302 339.2423 359.9579 349.2520
## 305 339.2423 359.9579 349.2520
## 307 339.2423 359.9579 349.2520
## 309 339.2423 359.9579 349.2520
## 311 339.2423 359.9579 349.2520
## 313 339.2423 359.9579 349.2520
## 315 339.2423 359.9579 349.2520
## 317 339.2423 359.9579 349.2520
## 319 339.2423 359.9579 349.2520
## 321 339.2423 359.9579 349.2520
## 323 339.2423 359.9579 349.2520
## 325 339.2423 359.9579 349.2520
## 327 339.2423 359.9579 349.2520
## 329 339.2423 359.9579 349.2520
## 331 339.2423 359.9579 349.2520
## 333 339.2423 359.9579 349.2520
## 335 339.2423 359.9579 349.2520
## 337 339.2423 359.9579 349.2520
## 339 339.2423 359.9579 349.2520
## 341 339.2423 359.9579 349.2520
## 343 339.2423 359.9579 349.2520
## 345 339.2423 359.9579 349.2520
## 347 339.2423 359.9579 349.2520
## 349 339.2423 359.9579 349.2520
## 351 339.2423 359.9579 349.2520
## 353 339.2423 359.9579 349.2520
## 355 339.2423 359.9579 349.2520
## 357 339.2423 359.9579 349.2520
## 359 339.2423 359.9579 349.2520
## 361 339.2423 359.9579 349.2520
## 363 339.2423 359.9579 349.2520
## 365 339.2423 359.9579 349.2520
## 367 339.2423 359.9579 349.2520
## 369 339.2423 359.9579 349.2520
## 371 339.2423 359.9579 349.2520
## 373 339.2423 359.9579 349.2520
## 375 339.2423 359.9579 349.2520
## 377 339.2423 359.9579 349.2520
## 379 339.2423 359.9579 349.2520
## 381 339.2423 359.9579 349.2520
## 383 339.2423 359.9579 349.2520
## 385 339.2423 359.9579 349.2520
## 387 339.2423 359.9579 349.2520
## 389 339.2423 359.9579 349.2520
## 391 339.2423 359.9579 349.2520
## 393 339.2423 359.9579 349.2520
## 395 339.2423 359.9579 349.2520
## 397 339.2423 359.9579 349.2520
## 399 339.2423 359.9579 349.2520
## 401 339.2423 359.9579 349.2520
## 403 339.2423 359.9579 349.2520
### Make frequency plots for paper


# predictions are from this model
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ (1 | beeID) + treatment
##    Data: newdf3
## 
## REML criterion at convergence: 25129.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3673 -0.4740  0.1874  0.6837  2.3307 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  901.5   30.02   
##  Residual             2357.8   48.56   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)        346.642      5.230   66.28
## treatmentWeighted    2.610      2.064    1.26
## 
## Correlation of Fixed Effects:
##             (Intr)
## trtmntWghtd -0.214
anova(m1, update(m1, .~. - treatment)) #p-val and BIC for paper
## refitting model(s) with ML (instead of REML)
## Data: newdf3
## Models:
## update(m1, . ~ . - treatment): freq ~ (1 | beeID)
## m1: freq ~ (1 | beeID) + treatment
##                               Df   AIC   BIC logLik deviance  Chisq Chi Df
## update(m1, . ~ . - treatment)  3 25146 25163 -12570    25140              
## m1                             4 25146 25169 -12569    25138 1.5992      1
##                               Pr(>Chisq)
## update(m1, . ~ . - treatment)           
## m1                                 0.206
pframe$trt3 <- mapvalues(pframe$trt2, from = c("Sham first", "Weighted first"), 
                          to = c("Sham -> Weighted", "Weighted -> Sham"))


g0 <- ggplot(pframe, aes(x=visitNum_centered, y=predMean))+
     geom_line(aes(color = treatment))+ 
     labs(y = "Frequency (Hz)", x = "Sonication number\n (0 is when treatment switched)") + 
     geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) + 
  facet_wrap(~trt3) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  theme(legend.position = c(0.5, 0.90), 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="horizontal") + 
  ylim(c(330, 370))
g0

ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries.pdf"), width =5, height = 3.6)


g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
                aes(x = visitNum_centered, y = freq), position = position_jitter(height = 5), size = 0.1, pch = 16, alpha = 0.4) + 
  ylim(c(220, 470))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
g1 # note that the distribution of frequency is a bit bimodal -- this could be 

# a problem, so we ran the analysis with a freq cutoff at 270
# we found no significant differences in the models either way.
ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries_RawData.pdf"), width =5, height = 3.6)


g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
                aes(x = visitNum_centered, y = freq), position = position_jitter(height = 5), size = 0.1, pch = 16, alpha = 0.4) + 
  ylim(c(270, 480))
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
g1
## Warning: Removed 730 rows containing missing values (geom_point).

ggsave(filename = file.path(figDir, "WeightedSham_Freq_Timeseries_RawData_trimmed.pdf"), width =5, height = 3.6)
## Warning: Removed 738 rows containing missing values (geom_point).
# mean and 95% CI (without visit number)

g0 <- ggplot(pframe[pframe$visitNum_centered %in% c(-25, 25), ], aes(x=treatment, y=predMean))+
    geom_point(aes(color = treatment), alpha = 1) + 
     labs(y = "Frequency (Hz)", x = "Flower treatment") + 
     geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  theme(legend.position = "none", 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="horizontal") 
g0

ggsave(filename = file.path(figDir, "WeightedSham_Freq_pointWhiskers.pdf"), width =5, height = 3.6)




## add raw data
g1 <- g0 + geom_point(data = newdf3, aes(x = treatment, y = freq, color = treatment), position = position_jitter(width = 0.1, height = 6), alpha = 0.4, pch = 16, size = 0.2)
g1

ggsave(filename = file.path(figDir, "WeightedSham_Freq_pointWhiskers_rawData.pdf"), width =5, height = 3.6)




# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
g0 <- ggplot(pframe, aes(x=visitNum_centered, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "buzzNumber") + 
     geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi), alpha = 0.2) + 
  facet_wrap(~trt2)
g0

# plot effect to see if it agrees
plot(Effect(c("treatment"), m1))

Modeling amplitude with gamm

g01 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2) + trt2 + IT + hive + treatment, random =  ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g01$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)

summary(g01$gam) # Summary for paper 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + trt2 + IT + 
##     hive + treatment
## 
## Parametric coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.72977    0.58219   4.689 2.90e-06 ***
## trt2Weighted first -0.05297    0.09833  -0.539    0.590    
## IT                  0.21347    0.14365   1.486    0.137    
## hive4               0.09955    0.10391   0.958    0.338    
## treatmentWeighted  -0.32857    0.04921  -6.678 3.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                           edf Ref.df     F p-value  
## s(visitNum_centered):trt2Sham first     1.414  1.414 4.924  0.0175 *
## s(visitNum_centered):trt2Weighted first 2.800  2.800 3.683  0.0323 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0571   
## lmer.REML = 4557.5  Scale est. = 0.384     n = 2360
summary(g01$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4557.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9897 -0.5342  0.1125  0.6707  2.6828 
## 
## Random effects:
##  Groups   Name                                    Variance Std.Dev.
##  beeID    (Intercept)                             0.076179 0.27601 
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.376254 0.61340 
##  Xr       s(visitNum_centered):trt2Sham first     0.008701 0.09328 
##  Residual                                         0.384002 0.61968 
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                              Estimate Std. Error t value
## X(Intercept)                                 2.729772   0.582186   4.689
## Xtrt2Weighted first                         -0.052966   0.098331  -0.539
## XIT                                          0.213470   0.143651   1.486
## Xhive4                                       0.099547   0.103913   0.958
## XtreatmentWeighted                          -0.328574   0.049206  -6.678
## Xs(visitNum_centered):trt2Sham firstFx1      0.059922   0.036356   1.648
## Xs(visitNum_centered):trt2Weighted firstFx1  0.003639   0.188745   0.019
## 
## Correlation of Fixed Effects:
##             X(Int) Xtr2Wf XIT    Xhive4 XtrtmW X(N_):2Sf
## Xtrt2Wghtdf -0.238                                      
## XIT         -0.991  0.166                               
## Xhive4      -0.041 -0.011 -0.012                        
## XtrtmntWght -0.062  0.061  0.015  0.000                 
## X(N_):2SfF1  0.051 -0.029 -0.030 -0.003 -0.455          
## X(N_):2WfF1 -0.001  0.060  0.000  0.001  0.023 -0.011
dev.off()
## null device 
##           1
g011 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2)  + IT + hive + treatment, random =  ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g011$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g011$gam) 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + hive + 
##     treatment
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.65539    0.55974   4.744 2.22e-06 ***
## IT                 0.22614    0.14023   1.613    0.107    
## hive4              0.09890    0.10286   0.961    0.336    
## treatmentWeighted -0.32629    0.04905  -6.652 3.58e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                           edf Ref.df     F p-value  
## s(visitNum_centered):trt2Sham first     1.396  1.396 4.868  0.0185 *
## s(visitNum_centered):trt2Weighted first 2.824  2.824 3.711  0.0286 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0543   
## lmer.REML =   4555  Scale est. = 0.38399   n = 2360
summary(g011$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4555
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9904 -0.5310  0.1120  0.6721  2.6775 
## 
## Random effects:
##  Groups   Name                                    Variance Std.Dev.
##  beeID    (Intercept)                             0.074527 0.27300 
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.387226 0.62227 
##  Xr       s(visitNum_centered):trt2Sham first     0.008222 0.09068 
##  Residual                                         0.383989 0.61967 
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 2.65539    0.55974   4.744
## XIT                                          0.22614    0.14023   1.613
## Xhive4                                       0.09890    0.10286   0.961
## XtreatmentWeighted                          -0.32629    0.04905  -6.652
## Xs(visitNum_centered):trt2Sham firstFx1      0.05943    0.03583   1.659
## Xs(visitNum_centered):trt2Weighted firstFx1  0.01018    0.19078   0.053
## 
## Correlation of Fixed Effects:
##             X(Int) XIT    Xhive4 XtrtmW X(N_):2Sf
## XIT         -0.994                               
## Xhive4      -0.045 -0.011                        
## XtrtmntWght -0.049  0.005  0.001                 
## X(N_):2SfF1  0.046 -0.026 -0.003 -0.463          
## X(N_):2WfF1  0.014 -0.010  0.002  0.018 -0.008
dev.off()
## null device 
##           1
anova(g01$mer, g011$mer) # go with g011
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g011$mer: NULL
## g01$mer: NULL
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## g011$mer 10 4553.4 4611.1 -2266.7   4533.4                         
## g01$mer  11 4555.1 4618.5 -2266.5   4533.1 0.3497      1     0.5543
g02 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2)  + IT  + treatment, random =  ~(1|beeID), data = newdf3)
par(mfrow = c(2,3))
aab <- plot(g02$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
summary(g02$gam) # Summary for paper 
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum_centered, by = trt2) + IT + treatment
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.67960    0.55806   4.802 1.67e-06 ***
## IT                 0.22770    0.13995   1.627    0.104    
## treatmentWeighted -0.32686    0.04912  -6.655 3.52e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                           edf Ref.df     F p-value  
## s(visitNum_centered):trt2Sham first     1.412  1.412 4.868  0.0182 *
## s(visitNum_centered):trt2Weighted first 2.810  2.810 3.669  0.0310 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0534   
## lmer.REML = 4553.2  Scale est. = 0.384     n = 2360
summary(g02$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4553.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9949 -0.5328  0.1143  0.6739  2.6875 
## 
## Random effects:
##  Groups   Name                                    Variance Std.Dev.
##  beeID    (Intercept)                             0.074211 0.27242 
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.379482 0.61602 
##  Xr       s(visitNum_centered):trt2Sham first     0.008663 0.09308 
##  Residual                                         0.384001 0.61968 
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                              Estimate Std. Error t value
## X(Intercept)                                 2.679603   0.558065   4.802
## XIT                                          0.227704   0.139951   1.627
## XtreatmentWeighted                          -0.326857   0.049116  -6.655
## Xs(visitNum_centered):trt2Sham firstFx1      0.059463   0.036299   1.638
## Xs(visitNum_centered):trt2Weighted firstFx1  0.009651   0.189097   0.051
## 
## Correlation of Fixed Effects:
##             X(Int) XIT    XtrtmW X(N_):2Sf
## XIT         -0.995                        
## XtrtmntWght -0.049  0.005                 
## X(N_):2SfF1  0.046 -0.026 -0.455          
## X(N_):2WfF1  0.014 -0.010  0.019 -0.009
dev.off()
## null device 
##           1
anova(g011$mer, g02$mer) # g02 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g02$mer: NULL
## g011$mer: NULL
##          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## g02$mer   9 4552.4 4604.3 -2267.2   4534.4                         
## g011$mer 10 4553.4 4611.1 -2266.7   4533.4 0.9778      1     0.3228
g03 = gamm4(log(amp_acc) ~ s(visitNum_centered, by = trt2)  + treatment, random =  ~(1|beeID), data = newdf3)

anova(g02$mer, g03$mer) #keep g03
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g03$mer: NULL
## g02$mer: NULL
##         Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## g03$mer  8 4553.1 4599.2 -2268.5   4537.1                        
## g02$mer  9 4552.4 4604.3 -2267.2   4534.4 2.691      1     0.1009
summary(g03$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4553.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9839 -0.5287  0.1117  0.6716  2.6967 
## 
## Random effects:
##  Groups   Name                                    Variance Std.Dev.
##  beeID    (Intercept)                             0.077961 0.27921 
##  Xr.0     s(visitNum_centered):trt2Weighted first 0.397553 0.63052 
##  Xr       s(visitNum_centered):trt2Sham first     0.005922 0.07696 
##  Residual                                         0.384011 0.61969 
## Number of obs: 2360, groups:  beeID, 36; Xr.0, 8; Xr, 8
## 
## Fixed effects:
##                                             Estimate Std. Error t value
## X(Intercept)                                 3.58211    0.05438   65.87
## XtreatmentWeighted                          -0.32459    0.04854   -6.69
## Xs(visitNum_centered):trt2Sham firstFx1      0.06166    0.03320    1.86
## Xs(visitNum_centered):trt2Weighted firstFx1  0.01321    0.19301    0.07
## 
## Correlation of Fixed Effects:
##             X(Int) XtrtmW X(N_):2Sf
## XtrtmntWght -0.451                 
## X(N_):2SfF1  0.228 -0.510          
## X(N_):2WfF1  0.040  0.017 -0.009
g04 <- gamm4(log(amp_acc) ~ s(visitNum_centered)  + treatment, random =  ~(1|beeID), data = newdf3)

anova(g03$mer, g04$mer) #g04 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g04$mer: NULL
## g03$mer: NULL
##         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## g04$mer  6 4551.2 4585.8 -2269.6   4539.2                         
## g03$mer  8 4553.1 4599.2 -2268.5   4537.1 2.1425      2     0.3426
summary(g04$gam) # approx p-value for smooth
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum_centered) + treatment
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.56539    0.05047   70.64   <2e-16 ***
## treatmentWeighted -0.29794    0.02688  -11.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df     F p-value   
## s(visitNum_centered) 2.056  2.056 6.862 0.00106 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0378   
## lmer.REML = 4555.5  Scale est. = 0.38524   n = 2360
summary(g04$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4555.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0433 -0.5427  0.1059  0.6832  2.7605 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev.
##  beeID    (Intercept)          0.07754  0.2785  
##  Xr       s(visitNum_centered) 0.02027  0.1424  
##  Residual                      0.38524  0.6207  
## Number of obs: 2360, groups:  beeID, 36; Xr, 8
## 
## Fixed effects:
##                          Estimate Std. Error t value
## X(Intercept)              3.56539    0.05047   70.64
## XtreatmentWeighted       -0.29794    0.02688  -11.08
## Xs(visitNum_centered)Fx1  0.02517    0.04290    0.59
## 
## Correlation of Fixed Effects:
##             X(Int) XtrtmW
## XtrtmntWght -0.293       
## Xs(vstN_)F1  0.024 -0.053
g44 <- gamm4(log(amp_acc) ~ visitNum_centered  + treatment, random =  ~(1|beeID), data = newdf3)
g05 <- gamm4(log(amp_acc) ~ treatment, random =  ~(1|beeID), data = newdf3)

anova(g44$mer, g05$mer) #g44 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g05$mer: NULL
## g44$mer: NULL
##         Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## g05$mer  4 4560.8 4583.9 -2276.4   4552.8                            
## g44$mer  5 4549.2 4578.0 -2269.6   4539.2 13.63      1  0.0002226 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
g06 <- gamm4(log(amp_acc) ~ visitNum_centered , random =  ~(1|beeID), data = newdf3)
anova(g44$mer, g06$mer) #g044 better
## refitting model(s) with ML (instead of REML)
## Data: NULL
## Models:
## g06$mer: NULL
## g44$mer: NULL
##         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## g06$mer  4 4666.2 4689.3 -2329.1   4658.2                             
## g44$mer  5 4549.2 4578.0 -2269.6   4539.2 119.02      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# final model -- only treatment matters
ampMod <- g04
summary(ampMod$mer)
## Linear mixed model fit by REML ['lmerMod']
## 
## REML criterion at convergence: 4555.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0433 -0.5427  0.1059  0.6832  2.7605 
## 
## Random effects:
##  Groups   Name                 Variance Std.Dev.
##  beeID    (Intercept)          0.07754  0.2785  
##  Xr       s(visitNum_centered) 0.02027  0.1424  
##  Residual                      0.38524  0.6207  
## Number of obs: 2360, groups:  beeID, 36; Xr, 8
## 
## Fixed effects:
##                          Estimate Std. Error t value
## X(Intercept)              3.56539    0.05047   70.64
## XtreatmentWeighted       -0.29794    0.02688  -11.08
## Xs(visitNum_centered)Fx1  0.02517    0.04290    0.59
## 
## Correlation of Fixed Effects:
##             X(Int) XtrtmW
## XtrtmntWght -0.293       
## Xs(vstN_)F1  0.024 -0.053
summary(ampMod$gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## log(amp_acc) ~ s(visitNum_centered) + treatment
## 
## Parametric coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.56539    0.05047   70.64   <2e-16 ***
## treatmentWeighted -0.29794    0.02688  -11.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df     F p-value   
## s(visitNum_centered) 2.056  2.056 6.862 0.00106 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0378   
## lmer.REML = 4555.5  Scale est. = 0.38524   n = 2360
plot(ampMod$mer)
par(mfrow = c(2,1))
plot(ampMod$gam, all.terms = TRUE, rug = FALSE, shade = TRUE)
dev.off()
## null device 
##           1

Modeling amplitude with lmer

# log transform for acceleration helps model fit better
# start with large model (including interaction)
ma0 = lmer(log(amp_acc) ~ visitNum_centered*trt2 +  hive +IT * treatment+ (1|beeID), data = newdf3, REML = FALSE)
summary(ma0)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## log(amp_acc) ~ visitNum_centered * trt2 + hive + IT * treatment +  
##     (1 | beeID)
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4541.7   4599.3  -2260.8   4521.7     2350 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0318 -0.5397  0.1116  0.6862  2.6657 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.06647  0.2578  
##  Residual             0.38294  0.6188  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                                        Estimate Std. Error t value
## (Intercept)                           3.3181167  0.5736148   5.785
## visitNum_centered                     0.0017926  0.0010587   1.693
## trt2Weighted first                   -0.0608425  0.0918382  -0.662
## hive4                                 0.0948197  0.0975308   0.972
## IT                                    0.0589548  0.1420758   0.415
## treatmentWeighted                    -1.4174293  0.3083754  -4.596
## visitNum_centered:trt2Weighted first -0.0004062  0.0020914  -0.194
## IT:treatmentWeighted                  0.2808255  0.0780233   3.599
## 
## Correlation of Fixed Effects:
##             (Intr) vstNm_ trt2Wf hive4  IT     trtmnW vN_:2f
## vstNm_cntrd -0.015                                          
## trt2Wghtdfr -0.229  0.040                                   
## hive4       -0.038 -0.014 -0.013                            
## IT          -0.992  0.038  0.162 -0.014                     
## trtmntWghtd -0.312  0.142  0.014 -0.009  0.316              
## vstNm_c:2Wf -0.022 -0.805 -0.027  0.005 -0.005 -0.002       
## IT:trtmntWg  0.307 -0.241 -0.016  0.010 -0.317 -0.990  0.115
ma1 = update(ma0, .~. - visitNum_centered:trt2)

BIC(ma0, ma1) # interaction, ma1, is better
##     df      BIC
## ma0 10 4599.347
## ma1  9 4591.618
# pval for interaction (in case paper needs it)
anova(ma0, ma1)
## Data: newdf3
## Models:
## ma1: log(amp_acc) ~ visitNum_centered + trt2 + hive + IT + treatment + 
## ma1:     (1 | beeID) + IT:treatment
## ma0: log(amp_acc) ~ visitNum_centered * trt2 + hive + IT * treatment + 
## ma0:     (1 | beeID)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ma1  9 4539.7 4591.6 -2260.9   4521.7                         
## ma0 10 4541.7 4599.3 -2260.8   4521.7 0.0377      1      0.846
summary(ma1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## log(amp_acc) ~ visitNum_centered + trt2 + hive + IT + treatment +  
##     (1 | beeID) + IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4539.7   4591.6  -2260.9   4521.7     2351 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0255 -0.5401  0.1127  0.6850  2.6689 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.06652  0.2579  
##  Residual             0.38294  0.6188  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)           3.3156701  0.5736617   5.780
## visitNum_centered     0.0016270  0.0006282   2.590
## trt2Weighted first   -0.0613216  0.0918373  -0.668
## hive4                 0.0949225  0.0975640   0.973
## IT                    0.0588142  0.1421196   0.414
## treatmentWeighted    -1.4175699  0.3083761  -4.597
## IT:treatmentWeighted  0.2825651  0.0775090   3.646
## 
## Correlation of Fixed Effects:
##             (Intr) vstNm_ trt2Wf hive4  IT     trtmnW
## vstNm_cntrd -0.055                                   
## trt2Wghtdfr -0.229  0.031                            
## hive4       -0.038 -0.016 -0.013                     
## IT          -0.993  0.056  0.162 -0.014              
## trtmntWghtd -0.312  0.236  0.014 -0.009  0.316       
## IT:trtmntWg  0.312 -0.253 -0.013  0.009 -0.318 -0.996
ma2 <- update(ma1, .~. - trt2)
BIC(ma1, ma2) # keep ma2
##     df      BIC
## ma1  9 4591.618
## ma2  8 4584.295
summary(ma2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(amp_acc) ~ visitNum_centered + hive + IT + treatment + (1 |  
##     beeID) + IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4538.2   4584.3  -2261.1   4522.2     2352 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0293 -0.5416  0.1160  0.6854  2.6651 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.06752  0.2598  
##  Residual             0.38294  0.6188  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           3.227873   0.561839   5.745
## visitNum_centered     0.001640   0.000628   2.611
## hive4                 0.094133   0.098224   0.958
## IT                    0.074144   0.141118   0.525
## treatmentWeighted    -1.414972   0.308358  -4.589
## IT:treatmentWeighted  0.281940   0.077505   3.638
## 
## Correlation of Fixed Effects:
##             (Intr) vstNm_ hive4  IT     trtmnW
## vstNm_cntrd -0.049                            
## hive4       -0.042 -0.016                     
## IT          -0.995  0.052 -0.012              
## trtmntWghtd -0.315  0.235 -0.008  0.316       
## IT:trtmntWg  0.316 -0.252  0.009 -0.318 -0.996
ma3 <- update(ma2, .~. - hive)
anova(ma2, ma3) # keep ma3
## Data: newdf3
## Models:
## ma3: log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) + 
## ma3:     IT:treatment
## ma2: log(amp_acc) ~ visitNum_centered + hive + IT + treatment + (1 | 
## ma2:     beeID) + IT:treatment
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## ma3  7 4537.1 4577.4 -2261.5   4523.1                         
## ma2  8 4538.2 4584.3 -2261.1   4522.2 0.9083      1     0.3406
summary(ma3)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) +  
##     IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4537.1   4577.4  -2261.5   4523.1     2353 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0319 -0.5419  0.1149  0.6887  2.6771 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.06929  0.2632  
##  Residual             0.38295  0.6188  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           3.250421   0.567438   5.728
## visitNum_centered     0.001648   0.000628   2.625
## IT                    0.075722   0.142637   0.531
## treatmentWeighted    -1.412882   0.308373  -4.582
## IT:treatmentWeighted  0.281383   0.077509   3.630
## 
## Correlation of Fixed Effects:
##             (Intr) vstNm_ IT     trtmnW
## vstNm_cntrd -0.049                     
## IT          -0.996  0.051              
## trtmntWghtd -0.312  0.235  0.312       
## IT:trtmntWg  0.313 -0.252 -0.315 -0.996
ma3.1 <- update(ma3, .~. - visitNum_centered)
anova(ma3, ma3.1) # keep 3.1
## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## ma3: log(amp_acc) ~ visitNum_centered + IT + treatment + (1 | beeID) + 
## ma3:     IT:treatment
##       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)   
## ma3.1  6 4541.9 4576.5 -2265.0   4529.9                           
## ma3    7 4537.1 4577.4 -2261.5   4523.1 6.877      1   0.008731 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ma3.1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4541.9   4576.5  -2265.0   4529.9     2354 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0202 -0.5393  0.1093  0.6834  2.6666 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.07006  0.2647  
##  Residual             0.38402  0.6197  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           3.32371    0.56952   5.836
## IT                    0.05660    0.14315   0.395
## treatmentWeighted    -1.60343    0.30014  -5.342
## IT:treatmentWeighted  0.33276    0.07511   4.431
## 
## Correlation of Fixed Effects:
##             (Intr) IT     trtmnW
## IT          -0.996              
## trtmntWghtd -0.309  0.308       
## IT:trtmntWg  0.310 -0.312 -0.996
ma4 <- update(ma3.1, .~. - visitNum_centered )
anova(ma3.1, ma4) #keep ma3.1
## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## ma4: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
##       Df    AIC    BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ma3.1  6 4541.9 4576.5  -2265   4529.9                        
## ma4    6 4541.9 4576.5  -2265   4529.9     0      0          1
summary(ma3.1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4541.9   4576.5  -2265.0   4529.9     2354 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0202 -0.5393  0.1093  0.6834  2.6666 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.07006  0.2647  
##  Residual             0.38402  0.6197  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           3.32371    0.56952   5.836
## IT                    0.05660    0.14315   0.395
## treatmentWeighted    -1.60343    0.30014  -5.342
## IT:treatmentWeighted  0.33276    0.07511   4.431
## 
## Correlation of Fixed Effects:
##             (Intr) IT     trtmnW
## IT          -0.996              
## trtmntWghtd -0.309  0.308       
## IT:trtmntWg  0.310 -0.312 -0.996
m3.2 <- update(ma3.1, .~ . + I(scale(visitNum_centered)) + I(scale(visitNum_centered^2)) + (1|beeID))
summary(m3.2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## log(amp_acc) ~ IT + treatment + (1 | beeID) + I(scale(visitNum_centered)) +  
##     I(scale(visitNum_centered^2)) + IT:treatment
##    Data: newdf3
## 
##      AIC      BIC   logLik deviance df.resid 
##   4539.1   4585.2  -2261.5   4523.1     2352 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0308 -0.5413  0.1154  0.6887  2.6758 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.06929  0.2632  
##  Residual             0.38294  0.6188  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                               Estimate Std. Error t value
## (Intercept)                    3.26035    0.56770   5.743
## IT                             0.07536    0.14273   0.528
## treatmentWeighted             -1.41125    0.30921  -4.564
## I(scale(visitNum_centered))    0.03784    0.01705   2.219
## I(scale(visitNum_centered^2))  0.00129    0.01797   0.072
## IT:treatmentWeighted           0.28096    0.07774   3.614
## 
## Correlation of Fixed Effects:
##             (Intr) IT     trtmnW I((N_) I((N_^
## IT          -0.996                            
## trtmntWghtd -0.307  0.309                     
## I(scl(vN_)) -0.057  0.062  0.164              
## I(sc(N_^2))  0.038 -0.035  0.074 -0.511       
## IT:trtmntWg  0.308 -0.311 -0.996 -0.177 -0.076
anova(ma3.1, m3.2) # keep ma3.1 (according to BIC)
## Data: newdf3
## Models:
## ma3.1: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
## m3.2: log(amp_acc) ~ IT + treatment + (1 | beeID) + I(scale(visitNum_centered)) + 
## m3.2:     I(scale(visitNum_centered^2)) + IT:treatment
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## ma3.1  6 4541.9 4576.5 -2265.0   4529.9                           
## m3.2   8 4539.1 4585.2 -2261.5   4523.1 6.8821      2    0.03203 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
finmod <- update(ma3.1, .~., REML = TRUE)
summary(finmod)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ IT + treatment + (1 | beeID) + IT:treatment
##    Data: newdf3
## 
## REML criterion at convergence: 4545.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0179 -0.5395  0.1071  0.6826  2.6697 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.07463  0.2732  
##  Residual             0.38434  0.6200  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           3.32401    0.58488   5.683
## IT                    0.05653    0.14700   0.385
## treatmentWeighted    -1.60422    0.30031  -5.342
## IT:treatmentWeighted  0.33294    0.07515   4.430
## 
## Correlation of Fixed Effects:
##             (Intr) IT     trtmnW
## IT          -0.996              
## trtmntWghtd -0.301  0.300       
## IT:trtmntWg  0.302 -0.303 -0.996
# diagnostics
plot(finmod)

qqnorm(ranef(finmod)$beeID[[1]])
qqline(ranef(finmod)$beeID[[1]])

sjp.lmer(finmod) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(finmod, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

Bootstrap CI’s for amplitude for figure for paper

# set number of bootstrap replicates for models
nbootSims2 = 1000

# using hive 3, since it's the one with the most data
pframe <- expand.grid(trt2 = levels(droplevels(newdf3$trt2)),
                     IT = ITmean,
                     visitNum_centered = -50:50, 
                     hive = factor(3, levels = levels(newdf3$hive)),  
                     beeID = 99999, 
                     treatment = levels(droplevels(newdf3$treatment)))
pframe$amp_acc <- 0


pframe2 <- pframe[(pframe$trt2 == "Sham first" &
                    pframe$treatment == "Sham" & 
                    pframe$visitNum_centered < 0) |
                    (pframe$trt2 == "Sham first" &
                    pframe$treatment == "Weighted" & 
                    pframe$visitNum_centered > 0) |
                    (pframe$trt2 == "Weighted first" &
                    pframe$treatment == "Weighted" & 
                    pframe$visitNum_centered < 0) |
                    (pframe$trt2 == "Weighted first" &
                    pframe$treatment == "Sham" & 
                    pframe$visitNum_centered > 0) 
                    , ]
# exponentiate to put on original scale
pp <- exp(predict(finmod, newdata = pframe2, re.form=NA, type = 'response')) # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(finmod, FUN=function(x) predict(x, pframe2, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe2$blo<- exp(bb2_se[1,]) # exponentiate to put on original scale
pframe2$bhi<- exp(bb2_se[2,])
pframe2$predMean <- pp
pframe2
##               trt2       IT visitNum_centered hive beeID treatment amp_acc
## 1       Sham first 3.965833               -50    3 99999      Sham       0
## 3       Sham first 3.965833               -49    3 99999      Sham       0
## 5       Sham first 3.965833               -48    3 99999      Sham       0
## 7       Sham first 3.965833               -47    3 99999      Sham       0
## 9       Sham first 3.965833               -46    3 99999      Sham       0
## 11      Sham first 3.965833               -45    3 99999      Sham       0
## 13      Sham first 3.965833               -44    3 99999      Sham       0
## 15      Sham first 3.965833               -43    3 99999      Sham       0
## 17      Sham first 3.965833               -42    3 99999      Sham       0
## 19      Sham first 3.965833               -41    3 99999      Sham       0
## 21      Sham first 3.965833               -40    3 99999      Sham       0
## 23      Sham first 3.965833               -39    3 99999      Sham       0
## 25      Sham first 3.965833               -38    3 99999      Sham       0
## 27      Sham first 3.965833               -37    3 99999      Sham       0
## 29      Sham first 3.965833               -36    3 99999      Sham       0
## 31      Sham first 3.965833               -35    3 99999      Sham       0
## 33      Sham first 3.965833               -34    3 99999      Sham       0
## 35      Sham first 3.965833               -33    3 99999      Sham       0
## 37      Sham first 3.965833               -32    3 99999      Sham       0
## 39      Sham first 3.965833               -31    3 99999      Sham       0
## 41      Sham first 3.965833               -30    3 99999      Sham       0
## 43      Sham first 3.965833               -29    3 99999      Sham       0
## 45      Sham first 3.965833               -28    3 99999      Sham       0
## 47      Sham first 3.965833               -27    3 99999      Sham       0
## 49      Sham first 3.965833               -26    3 99999      Sham       0
## 51      Sham first 3.965833               -25    3 99999      Sham       0
## 53      Sham first 3.965833               -24    3 99999      Sham       0
## 55      Sham first 3.965833               -23    3 99999      Sham       0
## 57      Sham first 3.965833               -22    3 99999      Sham       0
## 59      Sham first 3.965833               -21    3 99999      Sham       0
## 61      Sham first 3.965833               -20    3 99999      Sham       0
## 63      Sham first 3.965833               -19    3 99999      Sham       0
## 65      Sham first 3.965833               -18    3 99999      Sham       0
## 67      Sham first 3.965833               -17    3 99999      Sham       0
## 69      Sham first 3.965833               -16    3 99999      Sham       0
## 71      Sham first 3.965833               -15    3 99999      Sham       0
## 73      Sham first 3.965833               -14    3 99999      Sham       0
## 75      Sham first 3.965833               -13    3 99999      Sham       0
## 77      Sham first 3.965833               -12    3 99999      Sham       0
## 79      Sham first 3.965833               -11    3 99999      Sham       0
## 81      Sham first 3.965833               -10    3 99999      Sham       0
## 83      Sham first 3.965833                -9    3 99999      Sham       0
## 85      Sham first 3.965833                -8    3 99999      Sham       0
## 87      Sham first 3.965833                -7    3 99999      Sham       0
## 89      Sham first 3.965833                -6    3 99999      Sham       0
## 91      Sham first 3.965833                -5    3 99999      Sham       0
## 93      Sham first 3.965833                -4    3 99999      Sham       0
## 95      Sham first 3.965833                -3    3 99999      Sham       0
## 97      Sham first 3.965833                -2    3 99999      Sham       0
## 99      Sham first 3.965833                -1    3 99999      Sham       0
## 104 Weighted first 3.965833                 1    3 99999      Sham       0
## 106 Weighted first 3.965833                 2    3 99999      Sham       0
## 108 Weighted first 3.965833                 3    3 99999      Sham       0
## 110 Weighted first 3.965833                 4    3 99999      Sham       0
## 112 Weighted first 3.965833                 5    3 99999      Sham       0
## 114 Weighted first 3.965833                 6    3 99999      Sham       0
## 116 Weighted first 3.965833                 7    3 99999      Sham       0
## 118 Weighted first 3.965833                 8    3 99999      Sham       0
## 120 Weighted first 3.965833                 9    3 99999      Sham       0
## 122 Weighted first 3.965833                10    3 99999      Sham       0
## 124 Weighted first 3.965833                11    3 99999      Sham       0
## 126 Weighted first 3.965833                12    3 99999      Sham       0
## 128 Weighted first 3.965833                13    3 99999      Sham       0
## 130 Weighted first 3.965833                14    3 99999      Sham       0
## 132 Weighted first 3.965833                15    3 99999      Sham       0
## 134 Weighted first 3.965833                16    3 99999      Sham       0
## 136 Weighted first 3.965833                17    3 99999      Sham       0
## 138 Weighted first 3.965833                18    3 99999      Sham       0
## 140 Weighted first 3.965833                19    3 99999      Sham       0
## 142 Weighted first 3.965833                20    3 99999      Sham       0
## 144 Weighted first 3.965833                21    3 99999      Sham       0
## 146 Weighted first 3.965833                22    3 99999      Sham       0
## 148 Weighted first 3.965833                23    3 99999      Sham       0
## 150 Weighted first 3.965833                24    3 99999      Sham       0
## 152 Weighted first 3.965833                25    3 99999      Sham       0
## 154 Weighted first 3.965833                26    3 99999      Sham       0
## 156 Weighted first 3.965833                27    3 99999      Sham       0
## 158 Weighted first 3.965833                28    3 99999      Sham       0
## 160 Weighted first 3.965833                29    3 99999      Sham       0
## 162 Weighted first 3.965833                30    3 99999      Sham       0
## 164 Weighted first 3.965833                31    3 99999      Sham       0
## 166 Weighted first 3.965833                32    3 99999      Sham       0
## 168 Weighted first 3.965833                33    3 99999      Sham       0
## 170 Weighted first 3.965833                34    3 99999      Sham       0
## 172 Weighted first 3.965833                35    3 99999      Sham       0
## 174 Weighted first 3.965833                36    3 99999      Sham       0
## 176 Weighted first 3.965833                37    3 99999      Sham       0
## 178 Weighted first 3.965833                38    3 99999      Sham       0
## 180 Weighted first 3.965833                39    3 99999      Sham       0
## 182 Weighted first 3.965833                40    3 99999      Sham       0
## 184 Weighted first 3.965833                41    3 99999      Sham       0
## 186 Weighted first 3.965833                42    3 99999      Sham       0
## 188 Weighted first 3.965833                43    3 99999      Sham       0
## 190 Weighted first 3.965833                44    3 99999      Sham       0
## 192 Weighted first 3.965833                45    3 99999      Sham       0
## 194 Weighted first 3.965833                46    3 99999      Sham       0
## 196 Weighted first 3.965833                47    3 99999      Sham       0
## 198 Weighted first 3.965833                48    3 99999      Sham       0
## 200 Weighted first 3.965833                49    3 99999      Sham       0
## 202 Weighted first 3.965833                50    3 99999      Sham       0
## 204 Weighted first 3.965833               -50    3 99999  Weighted       0
## 206 Weighted first 3.965833               -49    3 99999  Weighted       0
## 208 Weighted first 3.965833               -48    3 99999  Weighted       0
## 210 Weighted first 3.965833               -47    3 99999  Weighted       0
## 212 Weighted first 3.965833               -46    3 99999  Weighted       0
## 214 Weighted first 3.965833               -45    3 99999  Weighted       0
## 216 Weighted first 3.965833               -44    3 99999  Weighted       0
## 218 Weighted first 3.965833               -43    3 99999  Weighted       0
## 220 Weighted first 3.965833               -42    3 99999  Weighted       0
## 222 Weighted first 3.965833               -41    3 99999  Weighted       0
## 224 Weighted first 3.965833               -40    3 99999  Weighted       0
## 226 Weighted first 3.965833               -39    3 99999  Weighted       0
## 228 Weighted first 3.965833               -38    3 99999  Weighted       0
## 230 Weighted first 3.965833               -37    3 99999  Weighted       0
## 232 Weighted first 3.965833               -36    3 99999  Weighted       0
## 234 Weighted first 3.965833               -35    3 99999  Weighted       0
## 236 Weighted first 3.965833               -34    3 99999  Weighted       0
## 238 Weighted first 3.965833               -33    3 99999  Weighted       0
## 240 Weighted first 3.965833               -32    3 99999  Weighted       0
## 242 Weighted first 3.965833               -31    3 99999  Weighted       0
## 244 Weighted first 3.965833               -30    3 99999  Weighted       0
## 246 Weighted first 3.965833               -29    3 99999  Weighted       0
## 248 Weighted first 3.965833               -28    3 99999  Weighted       0
## 250 Weighted first 3.965833               -27    3 99999  Weighted       0
## 252 Weighted first 3.965833               -26    3 99999  Weighted       0
## 254 Weighted first 3.965833               -25    3 99999  Weighted       0
## 256 Weighted first 3.965833               -24    3 99999  Weighted       0
## 258 Weighted first 3.965833               -23    3 99999  Weighted       0
## 260 Weighted first 3.965833               -22    3 99999  Weighted       0
## 262 Weighted first 3.965833               -21    3 99999  Weighted       0
## 264 Weighted first 3.965833               -20    3 99999  Weighted       0
## 266 Weighted first 3.965833               -19    3 99999  Weighted       0
## 268 Weighted first 3.965833               -18    3 99999  Weighted       0
## 270 Weighted first 3.965833               -17    3 99999  Weighted       0
## 272 Weighted first 3.965833               -16    3 99999  Weighted       0
## 274 Weighted first 3.965833               -15    3 99999  Weighted       0
## 276 Weighted first 3.965833               -14    3 99999  Weighted       0
## 278 Weighted first 3.965833               -13    3 99999  Weighted       0
## 280 Weighted first 3.965833               -12    3 99999  Weighted       0
## 282 Weighted first 3.965833               -11    3 99999  Weighted       0
## 284 Weighted first 3.965833               -10    3 99999  Weighted       0
## 286 Weighted first 3.965833                -9    3 99999  Weighted       0
## 288 Weighted first 3.965833                -8    3 99999  Weighted       0
## 290 Weighted first 3.965833                -7    3 99999  Weighted       0
## 292 Weighted first 3.965833                -6    3 99999  Weighted       0
## 294 Weighted first 3.965833                -5    3 99999  Weighted       0
## 296 Weighted first 3.965833                -4    3 99999  Weighted       0
## 298 Weighted first 3.965833                -3    3 99999  Weighted       0
## 300 Weighted first 3.965833                -2    3 99999  Weighted       0
## 302 Weighted first 3.965833                -1    3 99999  Weighted       0
## 305     Sham first 3.965833                 1    3 99999  Weighted       0
## 307     Sham first 3.965833                 2    3 99999  Weighted       0
## 309     Sham first 3.965833                 3    3 99999  Weighted       0
## 311     Sham first 3.965833                 4    3 99999  Weighted       0
## 313     Sham first 3.965833                 5    3 99999  Weighted       0
## 315     Sham first 3.965833                 6    3 99999  Weighted       0
## 317     Sham first 3.965833                 7    3 99999  Weighted       0
## 319     Sham first 3.965833                 8    3 99999  Weighted       0
## 321     Sham first 3.965833                 9    3 99999  Weighted       0
## 323     Sham first 3.965833                10    3 99999  Weighted       0
## 325     Sham first 3.965833                11    3 99999  Weighted       0
## 327     Sham first 3.965833                12    3 99999  Weighted       0
## 329     Sham first 3.965833                13    3 99999  Weighted       0
## 331     Sham first 3.965833                14    3 99999  Weighted       0
## 333     Sham first 3.965833                15    3 99999  Weighted       0
## 335     Sham first 3.965833                16    3 99999  Weighted       0
## 337     Sham first 3.965833                17    3 99999  Weighted       0
## 339     Sham first 3.965833                18    3 99999  Weighted       0
## 341     Sham first 3.965833                19    3 99999  Weighted       0
## 343     Sham first 3.965833                20    3 99999  Weighted       0
## 345     Sham first 3.965833                21    3 99999  Weighted       0
## 347     Sham first 3.965833                22    3 99999  Weighted       0
## 349     Sham first 3.965833                23    3 99999  Weighted       0
## 351     Sham first 3.965833                24    3 99999  Weighted       0
## 353     Sham first 3.965833                25    3 99999  Weighted       0
## 355     Sham first 3.965833                26    3 99999  Weighted       0
## 357     Sham first 3.965833                27    3 99999  Weighted       0
## 359     Sham first 3.965833                28    3 99999  Weighted       0
## 361     Sham first 3.965833                29    3 99999  Weighted       0
## 363     Sham first 3.965833                30    3 99999  Weighted       0
## 365     Sham first 3.965833                31    3 99999  Weighted       0
## 367     Sham first 3.965833                32    3 99999  Weighted       0
## 369     Sham first 3.965833                33    3 99999  Weighted       0
## 371     Sham first 3.965833                34    3 99999  Weighted       0
## 373     Sham first 3.965833                35    3 99999  Weighted       0
## 375     Sham first 3.965833                36    3 99999  Weighted       0
## 377     Sham first 3.965833                37    3 99999  Weighted       0
## 379     Sham first 3.965833                38    3 99999  Weighted       0
## 381     Sham first 3.965833                39    3 99999  Weighted       0
## 383     Sham first 3.965833                40    3 99999  Weighted       0
## 385     Sham first 3.965833                41    3 99999  Weighted       0
## 387     Sham first 3.965833                42    3 99999  Weighted       0
## 389     Sham first 3.965833                43    3 99999  Weighted       0
## 391     Sham first 3.965833                44    3 99999  Weighted       0
## 393     Sham first 3.965833                45    3 99999  Weighted       0
## 395     Sham first 3.965833                46    3 99999  Weighted       0
## 397     Sham first 3.965833                47    3 99999  Weighted       0
## 399     Sham first 3.965833                48    3 99999  Weighted       0
## 401     Sham first 3.965833                49    3 99999  Weighted       0
## 403     Sham first 3.965833                50    3 99999  Weighted       0
##          blo      bhi predMean
## 1   27.88884 35.12006 31.19904
## 3   27.97216 35.17425 31.26938
## 5   28.04596 35.23476 31.33989
## 7   28.12669 35.29929 31.41055
## 9   28.20839 35.33684 31.48137
## 11  28.29033 35.39071 31.55236
## 13  28.37251 35.44539 31.62350
## 15  28.45503 35.50759 31.69480
## 17  28.53821 35.57252 31.76626
## 19  28.59370 35.63709 31.83789
## 21  28.68210 35.70151 31.90967
## 23  28.75050 35.76053 31.98162
## 25  28.81907 35.81630 32.05373
## 27  28.88780 35.86353 32.12601
## 29  28.95670 35.90168 32.19844
## 31  29.02589 35.93960 32.27104
## 33  29.11343 35.99425 32.34380
## 35  29.20606 36.06262 32.41673
## 37  29.29943 36.10986 32.48982
## 39  29.39076 36.16208 32.56308
## 41  29.48507 36.21724 32.63650
## 43  29.56848 36.27248 32.71008
## 45  29.64726 36.32781 32.78384
## 47  29.72670 36.38323 32.85776
## 49  29.80635 36.45501 32.93184
## 51  29.87136 36.54930 33.00609
## 53  29.91921 36.64020 33.08051
## 55  29.99823 36.72050 33.15510
## 57  30.09295 36.77759 33.22986
## 59  30.18797 36.83477 33.30478
## 61  30.28329 36.89200 33.37988
## 63  30.35511 36.99573 33.45514
## 65  30.40840 37.11562 33.53057
## 67  30.46235 37.17469 33.60617
## 69  30.52589 37.20508 33.68195
## 71  30.56766 37.29571 33.75789
## 73  30.61995 37.37031 33.83401
## 75  30.67267 37.41339 33.91029
## 77  30.72607 37.47907 33.98675
## 79  30.80508 37.55789 34.06338
## 81  30.91853 37.62946 34.14019
## 83  31.03241 37.69277 34.21716
## 85  31.14670 37.77141 34.29432
## 87  31.23071 37.86450 34.37164
## 89  31.27387 37.95638 34.44914
## 91  31.36409 38.02615 34.52681
## 93  31.43114 38.11712 34.60466
## 95  31.51526 38.19815 34.68269
## 97  31.62292 38.27973 34.76089
## 99  31.71499 38.36068 34.83926
## 104 31.81108 38.52180 34.99655
## 106 31.83995 38.60391 35.07546
## 108 31.90001 38.70058 35.15454
## 110 31.96790 38.78534 35.23381
## 112 32.03549 38.85674 35.31325
## 114 32.10339 38.94356 35.39287
## 116 32.18359 39.04564 35.47267
## 118 32.23909 39.14751 35.55265
## 120 32.29222 39.22410 35.63282
## 122 32.33937 39.30174 35.71316
## 124 32.39874 39.39584 35.79368
## 126 32.47772 39.51021 35.87439
## 128 32.55723 39.60449 35.95527
## 130 32.63758 39.69793 36.03634
## 132 32.72030 39.79161 36.11760
## 134 32.79744 39.86463 36.19903
## 136 32.87484 39.95446 36.28065
## 138 32.93179 40.05022 36.36245
## 140 32.99830 40.13817 36.44444
## 142 33.08852 40.25472 36.52662
## 144 33.18291 40.37152 36.60897
## 146 33.23159 40.48940 36.69152
## 148 33.31052 40.59253 36.77425
## 150 33.37477 40.69400 36.85716
## 152 33.43175 40.79572 36.94027
## 154 33.49066 40.89770 37.02356
## 156 33.56110 40.99993 37.10704
## 158 33.62460 41.10234 37.19070
## 160 33.68271 41.23674 37.27456
## 162 33.73211 41.39638 37.35860
## 164 33.77496 41.52361 37.44284
## 166 33.83217 41.63195 37.52726
## 168 33.88948 41.73962 37.61187
## 170 33.94689 41.84628 37.69668
## 172 34.00439 41.95694 37.78167
## 174 34.06199 42.06651 37.86686
## 176 34.11969 42.17577 37.95224
## 178 34.17749 42.28531 38.03781
## 180 34.23539 42.40321 38.12358
## 182 34.28500 42.52645 38.20954
## 184 34.32780 42.65014 38.29569
## 186 34.37644 42.78285 38.38204
## 188 34.41201 42.89541 38.46858
## 190 34.45351 43.00827 38.55531
## 192 34.52197 43.15890 38.64225
## 194 34.56409 43.31342 38.72938
## 196 34.62097 43.45801 38.81670
## 198 34.66255 43.61892 38.90422
## 200 34.70441 43.77043 38.99194
## 202 34.74647 43.91361 39.07986
## 204 20.55672 26.07949 23.17974
## 206 20.61690 26.11577 23.23201
## 208 20.67727 26.15091 23.28439
## 210 20.73780 26.17571 23.33689
## 212 20.79077 26.20054 23.38951
## 214 20.85701 26.23645 23.44224
## 216 20.92717 26.27437 23.49510
## 218 20.99775 26.31200 23.54808
## 220 21.06808 26.34779 23.60117
## 222 21.13205 26.38790 23.65439
## 224 21.18395 26.42474 23.70772
## 226 21.23607 26.46435 23.76117
## 228 21.28905 26.50667 23.81475
## 230 21.36198 26.54996 23.86845
## 232 21.42087 26.59180 23.92226
## 234 21.47554 26.63316 23.97620
## 236 21.53767 26.67064 24.03026
## 238 21.60035 26.70818 24.08444
## 240 21.66341 26.76016 24.13875
## 242 21.72261 26.80296 24.19317
## 244 21.77216 26.84536 24.24772
## 246 21.83306 26.88782 24.30239
## 248 21.90161 26.93232 24.35719
## 250 21.94748 26.99432 24.41211
## 252 22.00056 27.04006 24.46715
## 254 22.06506 27.07071 24.52232
## 256 22.12992 27.10140 24.57761
## 258 22.19517 27.13248 24.63303
## 260 22.21556 27.16865 24.68857
## 262 22.30524 27.21512 24.74423
## 264 22.39241 27.26625 24.80003
## 266 22.45861 27.31336 24.85594
## 268 22.52501 27.36152 24.91199
## 270 22.57469 27.41286 24.96816
## 272 22.61547 27.46429 25.02445
## 274 22.67089 27.51562 25.08088
## 276 22.73343 27.56726 25.13743
## 278 22.77901 27.61896 25.19411
## 280 22.84112 27.64968 25.25091
## 282 22.88881 27.72290 25.30785
## 284 22.93692 27.77492 25.36491
## 286 22.98559 27.82703 25.42210
## 288 23.03671 27.87924 25.47942
## 290 23.08128 27.94271 25.53687
## 292 23.13520 27.99616 25.59445
## 294 23.18190 28.05627 25.65216
## 296 23.22688 28.12403 25.71000
## 298 23.27534 28.16709 25.76797
## 300 23.32393 28.22074 25.82607
## 302 23.37327 28.26137 25.88430
## 305 23.48898 28.43993 26.00115
## 307 23.53400 28.49247 26.05978
## 309 23.57911 28.54509 26.11854
## 311 23.62463 28.61448 26.17743
## 313 23.69228 28.69139 26.23645
## 315 23.76021 28.76851 26.29561
## 317 23.82834 28.81175 26.35490
## 319 23.89666 28.86456 26.41432
## 321 23.96518 28.93212 26.47388
## 323 24.01281 29.01442 26.53357
## 325 24.07782 29.06770 26.59339
## 327 24.12995 29.14005 26.65336
## 329 24.17227 29.21513 26.71345
## 331 24.22606 29.28164 26.77368
## 333 24.26679 29.35011 26.83405
## 335 24.30759 29.42384 26.89456
## 337 24.35299 29.48817 26.95520
## 339 24.41215 29.56041 27.01597
## 341 24.47146 29.63282 27.07689
## 343 24.51992 29.68743 27.13794
## 345 24.56797 29.75723 27.19913
## 347 24.61660 29.82718 27.26045
## 349 24.66580 29.90754 27.32192
## 351 24.73277 29.96894 27.38352
## 353 24.81203 30.04055 27.44526
## 355 24.85470 30.14112 27.50715
## 357 24.89734 30.21195 27.56917
## 359 24.94006 30.25096 27.63133
## 361 24.96859 30.32426 27.69363
## 363 25.00721 30.39774 27.75607
## 365 25.05565 30.47804 27.81866
## 367 25.10524 30.58094 27.88138
## 369 25.15427 30.66954 27.94424
## 371 25.20482 30.76536 28.00725
## 373 25.24056 30.83883 28.07040
## 375 25.27284 30.89840 28.13369
## 377 25.32440 30.95795 28.19713
## 379 25.37020 31.01827 28.26070
## 381 25.41419 31.09500 28.32442
## 383 25.45712 31.18639 28.38829
## 385 25.50723 31.26283 28.45230
## 387 25.56073 31.31911 28.51645
## 389 25.61107 31.42121 28.58074
## 391 25.65226 31.52449 28.64519
## 393 25.68123 31.62820 28.70977
## 395 25.73948 31.73151 28.77451
## 397 25.78224 31.80616 28.83939
## 399 25.82513 31.87966 28.90441
## 401 25.86815 31.97789 28.96958
## 403 25.91125 32.06483 29.03490
### Plot amplitude lmer with trial number

pframe2$trt3 <- mapvalues(pframe2$trt2, from = c("Sham first", "Weighted first"), 
                          to = c("Sham -> Weighted", "Weighted -> Sham"))

g0 <- ggplot(pframe2, aes(x=visitNum_centered, y=predMean))+
     geom_line(aes(color = treatment))+ 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication number\n (0 is when treatment switched)") + 
     geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) + 
  facet_wrap(~trt3) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  theme(legend.position = c(0.5, 0.89), 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="horizontal")
g0

ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_amp.pdf"), width =5, height = 3.6)
# I like this plot


g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
                aes(x = visitNum_centered, y = amp_acc), position = position_jitter(width = 0.5),  size = 0.1, pch = 16, alpha = 0.4) + 
  ylim(c(0, 165)) 
g1

ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_withRawData.pdf"), width =5, height = 3.6)
# I like this plot



g1 <- g0 + geom_point(data = newdf3[newdf3$visitNum_centered < 50,],
                aes(x = visitNum_centered, y = amp_acc), position = position_jitter(width = 0.5),  size = 0.1, pch = 16, alpha = 0.4) + 
  scale_y_log10(limits = c(2, 500)) + 
  annotation_logticks(sides = "l", color = 'grey', size = 0.2) 
g1

ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_withRawData_logScale.pdf"), width =5, height = 3.6)

# make figure that isn't faceted by order of treatment

g03 <- ggplot(pframe2, aes(x=visitNum_centered, y=predMean))+
     geom_line(aes(color = treatment))+ 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Sonication number\n (0 is when treatment switched)") + 
     geom_ribbon(aes(x = visitNum_centered, ymin = blo, ymax = bhi, fill= treatment), alpha = 0.2) + 
  facet_wrap(~treatment) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  scale_fill_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  theme(legend.position = c(0.75, 0.8), 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="vertical")
  
g03

ggsave(filename = file.path(figDir, "WeightedSham_Timeseries_trtFacet.pdf"), width =5, height = 3.6)
# I like this plot

Show amplitude plot without visit number

look at visit number 1 (right at the time of switching)

# mean and 95% CI (without visit number)

g0 <- ggplot(pframe2[pframe2$visitNum_centered %in% c(1), ], aes(x=treatment, y=predMean))+
    geom_point(aes(color = treatment), alpha = 1) + 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") + 
     geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, option = "A", begin =0, end = 0) + 
  theme(legend.position = "none", 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="horizontal") 
g0

ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers_BW.pdf"), width =5, height = 3.6)



g0 <- ggplot(pframe2[pframe2$visitNum_centered %in% c(1), ], aes(x=treatment, y=predMean))+
    geom_point(aes(color = treatment), alpha = 1) + 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") + 
     geom_errorbar(aes(x = treatment, ymin = blo, ymax = bhi, color= treatment), alpha = 1, width = 0.1) + 
  scale_color_viridis(name = "Flower treatment", discrete = TRUE, begin =0.3, end = 0.8) + 
  theme(legend.position = "none", 
        legend.background = element_rect(fill=alpha('gray95', 1)), 
        legend.direction="horizontal") 
g0

ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers.pdf"), width =5, height = 3.6)

g22 <- g0 + 
  annotate(geom="text", x=c(1,2), y=c(0, 0) + 40, label=c("a", "b"),
                color="black")
ggsave(filename = file.path(figDir, "WeightedSham_Amp_pointWhiskers_annot.pdf"), width =5, height = 3.6)

## add raw data
g1 <- g0 + geom_point(data = newdf3, aes(x = treatment, y = amp_acc, color = treatment), position = position_jitter(width = 0.1, height = 6), alpha = 0.4, pch = 16, size = 0.2)
g1

ggsave(filename = file.path(figDir, "WeightedSham_amp_pointWhiskers_rawData.pdf"), width =5, height = 3.6)

major break – below this line is the first method I used

this is almost identical to the code above

Modeling frequency with lmer

Use BIC to select model (decide what interactions and covariates to use)

Make sure that REML = FALSE when comparing BIC values

Start with the biggest model of interest, and then see what predictors can be removed

# interaction that may be important, based on domain knowledge
m0 = lmer(freq ~ treatment*IT+  hive + (1|beeID), data = new_df, REML = FALSE)

# no interaction
m1 = lmer(freq ~ treatment+IT+  hive + (1|beeID), data = new_df, REML = FALSE)
BIC(m0, m1) # stay with m1 -- no interaction
##    df      BIC
## m0  7 25187.77
## m1  6 25184.69
anova(m0, m1) # this p-value is very borderline -- we're being 
## Data: new_df
## Models:
## m1: freq ~ treatment + IT + hive + (1 | beeID)
## m0: freq ~ treatment * IT + hive + (1 | beeID)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m1  6 25150 25185 -12569    25138                           
## m0  7 25147 25188 -12567    25133 4.6848      1    0.03043 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# conservative to avoid overfitting by using BIC

m2 = lmer(freq ~ hive + IT+ (1|beeID), data = new_df, REML = FALSE)
BIC(m1, m2) # treatment doesn't make the model better
##    df      BIC
## m1  6 25184.69
## m2  5 25178.53
anova(m1, m2) # p-value for paper -- probably don't need to report
## Data: new_df
## Models:
## m2: freq ~ hive + IT + (1 | beeID)
## m1: freq ~ treatment + IT + hive + (1 | beeID)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2  5 25150 25178 -12570    25140                         
## m1  6 25150 25185 -12569    25138 1.6039      1     0.2054
# p-value for hive
m3 = lmer(freq ~ treatment + (1|beeID), data = new_df, REML = FALSE)
anova(m1, m3)
## Data: new_df
## Models:
## m3: freq ~ treatment + (1 | beeID)
## m1: freq ~ treatment + IT + hive + (1 | beeID)
##    Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  4 25146 25169 -12569    25138                         
## m1  6 25150 25185 -12569    25138 0.1592      2     0.9235
# diagnostics -- use REML = TRUE
m1 <- update(m1, .~., REML =TRUE)
summary(m1) # summary for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ treatment + IT + hive + (1 | beeID)
##    Data: new_df
## 
## REML criterion at convergence: 25115.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3676 -0.4684  0.1873  0.6825  2.3319 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept)  954.3   30.89   
##  Residual             2357.8   48.56   
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)        355.173     61.937   5.734
## treatmentWeighted    2.614      2.065   1.266
## IT                  -1.833     15.543  -0.118
## hive4               -4.129     11.398  -0.362
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnW IT    
## trtmntWghtd -0.002              
## IT          -0.995 -0.016       
## hive4       -0.047  0.000 -0.009
plot(m1)

qqnorm(ranef(m1)$beeID[[1]])
qqline(ranef(m1)$beeID[[1]])

sjp.lmer(m1, type = "re", sort = TRUE) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(m1, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

Bootstrap CI’s for figure for paper

# set number of bootstrap replicates for models
nbootSims = 1000

table(new_df$hive) # more trials from hive 3
## 
##    3    4 
## 1644  716
# using hive 3, since it's the one with the most data
# however, hive doesn't affect model anyway

# calculate an average IT for prediction
ITmean = mean(tapply(new_df$IT, INDEX = new_df$beeID, FUN = function(x) x[1] ))

pframe <- data.frame(treatment = levels(droplevels(new_df$treatment)),
                     IT = ITmean, 
                     hive = factor(3, levels = levels(new_df$hive)),  
                     beeID = 99999)
pframe$freq <- 0
pp <- predict(m1, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m1, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('treatment', "blo", "bhi", "predMean")]
pframe
##   treatment      blo      bhi predMean
## 1      Sham 335.3492 359.8953 347.9021
## 2  Weighted 337.7346 363.1810 350.5156

Make frequency plots for paper

# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
g0 <- ggplot(pframe, aes(x=treatment, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Flower treatment") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)
g0

ggsave(plot = g0, filename = file.path(figDir, "HeavyLight_PredsAndCI_freq.pdf"), width =4, height = 3)

Modeling amplitude with lmer

# log transform for acceleration helps model fit better
# start with large model (including interaction)
ma0 = lmer(log(amp_acc) ~ treatment * IT +  hive + (1|beeID), data = new_df, REML = FALSE)

ma1 = lmer(log(amp_acc) ~ treatment + IT +  hive + (1|beeID), data = new_df, REML = FALSE)

BIC(ma0, ma1) # interaction, ma0, is better
##     df      BIC
## ma0  7 4583.334
## ma1  6 4595.150
# pval for interaction (in case paper needs it)
anova(ma0, ma1)
## Data: new_df
## Models:
## ma1: log(amp_acc) ~ treatment + IT + hive + (1 | beeID)
## ma0: log(amp_acc) ~ treatment * IT + hive + (1 | beeID)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## ma1  6 4560.6 4595.1 -2274.3   4548.6                             
## ma0  7 4543.0 4583.3 -2264.5   4529.0 19.582      1  9.639e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ma2 <- update(ma0, .~. - hive)
BIC(ma0, ma2) # hive not important
##     df      BIC
## ma0  7 4583.334
## ma2  6 4576.547
#pval for hive
anova(ma0, ma2)
## Data: new_df
## Models:
## ma2: log(amp_acc) ~ treatment + IT + (1 | beeID) + treatment:IT
## ma0: log(amp_acc) ~ treatment * IT + hive + (1 | beeID)
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## ma2  6 4541.9 4576.5 -2265.0   4529.9                        
## ma0  7 4543.0 4583.3 -2264.5   4529.0 0.979      1     0.3224
# fit with reml = TRUE for summarizing
ma3 <- update(ma2,.~., REML = TRUE)
summary(ma3) # final model to report
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ treatment + IT + (1 | beeID) + treatment:IT
##    Data: new_df
## 
## REML criterion at convergence: 4545.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0179 -0.5395  0.1071  0.6826  2.6697 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeID    (Intercept) 0.07463  0.2732  
##  Residual             0.38434  0.6200  
## Number of obs: 2360, groups:  beeID, 36
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)           3.32401    0.58488   5.683
## treatmentWeighted    -1.60422    0.30031  -5.342
## IT                    0.05653    0.14700   0.385
## treatmentWeighted:IT  0.33294    0.07515   4.430
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnW IT    
## trtmntWghtd -0.301              
## IT          -0.996  0.300       
## trtmntWg:IT  0.302 -0.996 -0.303
# diagnostics
plot(ma3)

qqnorm(ranef(ma3)$beeID[[1]])
qqline(ranef(ma3)$beeID[[1]])

sjp.lmer(ma3) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(ma3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## `sjp.lmer()` will become deprecated in the future. Please use `plot_model()` instead.
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

Bootstrap CI’s for amplitude for figure for paper

# set number of bootstrap replicates for models
nbootSims2 = 1000

# using hive 3, since it's the one with the most data
pframe <- data.frame(treatment = levels(droplevels(new_df$treatment)), IT = ITmean, hive = factor(3, levels = levels(new_df$hive)),  beeID = 99999)
pframe$amp_acc <- 0
# exponentiate to put on original scale
pp <- exp(predict(ma3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(ma3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<- exp(bb2_se[1,]) # exponentiate to put on original scale
pframe$bhi<- exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('treatment', "blo", "bhi", "predMean")]
pframe
##   treatment      blo      bhi predMean
## 1      Sham 31.55900 37.97060 34.75029
## 2  Weighted 23.82485 28.84456 26.16290

Make amplitude plots for paper

# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# Holding IT = meanIT
print(ITmean)
## [1] 3.965833
ga0 <- ggplot(pframe, aes(x=treatment, y=predMean))+
     geom_point()+ 
     labs(y = expression ("Sonication amplitude "(m~s^{-2})), x = "Flower treatment") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1) +
  annotate(geom="text", x=c(1,2), y=c(0, 0) + 40, label=c("a", "b"),
                color="black") 
ga0

ggsave(plot = ga0, filename = file.path(figDir, "HeavyLight_PredsAndCI_amp.pdf"), width =4, height = 3)
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.4.0     viridisLite_0.2.0 gamm4_0.2-5      
##  [4] mgcv_1.8-20       nlme_3.1-131      effects_4.0-0    
##  [7] carData_3.0-0     plyr_1.8.4        multcomp_1.4-8   
## [10] TH.data_1.0-8     MASS_7.3-47       survival_2.41-3  
## [13] mvtnorm_1.0-6     sjPlot_2.4.0      lme4_1.1-14      
## [16] Matrix_1.2-11     reshape2_1.4.2    ggplot2_2.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] pbkrtest_0.4-7     RColorBrewer_1.1-2 rprojroot_1.2     
##  [4] tools_3.4.2        TMB_1.7.11         backports_1.1.1   
##  [7] R6_2.2.2           sjlabelled_1.0.5   DT_0.2            
## [10] lazyeval_0.2.1     colorspace_1.3-2   nnet_7.3-12       
## [13] gridExtra_2.3      tidyselect_0.2.3   mnormt_1.5-5      
## [16] compiler_3.4.2     sandwich_2.4-0     labeling_0.3      
## [19] scales_0.5.0       lmtest_0.9-35      psych_1.7.8       
## [22] blme_1.0-4         stringr_1.2.0      digest_0.6.12     
## [25] foreign_0.8-69     minqa_1.2.4        rmarkdown_1.8     
## [28] stringdist_0.9.4.6 pkgconfig_2.0.1    htmltools_0.3.6   
## [31] pwr_1.2-1          htmlwidgets_0.9    rlang_0.1.4       
## [34] shiny_1.0.5        bindr_0.1          zoo_1.8-0         
## [37] dplyr_0.7.4        magrittr_1.5       modeltools_0.2-21 
## [40] bayesplot_1.4.0    Rcpp_0.12.13       munsell_0.4.3     
## [43] abind_1.4-5        prediction_0.2.0   stringi_1.1.6     
## [46] yaml_2.1.14        merTools_0.3.0     snakecase_0.5.1   
## [49] grid_3.4.2         parallel_3.4.2     sjmisc_2.6.2      
## [52] forcats_0.2.0      lattice_0.20-35    ggeffects_0.2.2   
## [55] haven_1.1.0        splines_3.4.2      sjstats_0.12.0    
## [58] knitr_1.17         codetools_0.2-15   stats4_3.4.2      
## [61] glue_1.2.0         evaluate_0.10.1    modelr_0.1.1      
## [64] nloptr_1.0.4       httpuv_1.3.5       gtable_0.2.0      
## [67] purrr_0.2.4        tidyr_0.7.2        assertthat_0.2.0  
## [70] mime_0.5           coin_1.2-1         xtable_1.8-2      
## [73] broom_0.4.2        survey_3.32-1      coda_0.19-1       
## [76] tibble_1.3.4       arm_1.9-3          glmmTMB_0.1.4     
## [79] bindrcpp_0.2